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Within  the  last  decade  the  anticipation  of  some  experimentalists  for  the  ability  to 
ascertain  the  structure  of  biologically  important  chemical  systems  has  been  realized.  The 
elucidation  of  the  structure  of  the  photosynthetic  reaction  center  of  R.  Viridis,  including 
some  13,000  atomic  centers,  is  a  prime  example,  and  this  structure  has  led  to  a  deeper 
understanding  of  the  mechanisms  of  photosynthesis.  Along  with  this  increase  in  structural 
knowledge  a  concurrent  increase  in  the  knowledge,  of  the  observable  properties  of  such 
structures  has  also  occurred.  Of  particular  interest  to  this  work  are  the  magnetic  properties 
of  these  systems,  as  revealed  by  Electron  Spin  Resonance  (ESR)  spectroscopy.  The 
magnetic  spectra,  ESR,  of  such  systems  are  quite  complex  and  require  some  effort  to 
understand,  and  a  need  for  theoretical  aids  has  been  expressed.  In  many  cases  the  ESR 

vill 


signal  is  the  only  structural  signature  available. 

This  work  has  explored  several  quantum  chemical  methods  for  determining  the 
isotropic  hyperfine  splitting  (hfs)  in  such  systems  as  those  described  above.  Several 
semiempirical  Self-Consistent  Field  (SCF)  and  post-SCF  approaches  are  explored.  Two 
distinct  parameterizations  of  the  Intermediate  Neglect  of  Differential  Overlap  (INDO) 
approximation  are  examined  for  their  ability  to  produce  good  reference  functions.  Both 
Unrestricted  Hartree-Fock  (UHF)  and  Restricted  Open-Shell  Hartree-Fock  (ROHF)  SCF 
wavefunctions  are  investigated,  as  well  as  the  Configuration  Interaction  Single  replace- 
ment (CIS)  wavefunction.  The  results  of  applying  the  most  promising  method  to  several 
of  the  bacteria  chlorophyll  ions  are  reported  and  compared  to  those  observed  experi- 
mentally. 
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CHAPTER  1 
INTRODUCTION 

The  ESR  Experiment 

In  general,  magnetic  resonance  spectroscopy  observes  transitions  between  energy 
levels  which  depend  on  the  strength  of  a  magnetic  field  [1-3].  Experimental  techniques 
which  measure  the  interaction  between  an  external  magnetic  field  and  the  intrinsic  elec- 
tron spin  magnetic  moment  [4,  5]  are  usually  labeled  electron  spin  resonance  (ESR) 
spectroscopy  experiments.  Techniques  which  deal  with  interactions  between  orbital  mag- 
netic moments  and  an  external  magnetic  field  are  labeled  electron  paramagnetic  resonance 
(EPR)  specu-oscopy  experiments.  These  terms  are  sometime  used  interchangeably;  in  this 
work  we  shall  deal  with  the  ESR  problem  as  defined  above. 

G 

While  any  atomic  or  molecular  system  containing  one  or  more  unpaired  electrons 
will  in  general  possess  a  magnetic  moment  due  to  the  intrinsic  spin  of  electrons,  we  shall 
be  concerned  in  this  work  with  molecular  systems  containing  a  single  unpaired  electron. 
(And  as  we  shall  see  later,  really  only  with  those  magnetic  molecular  systems  that  are 
"large"  in  the  sense  of  ab-initio  quantum  chemical  theory.)  This  work  is  an  attempt  to 
develop  a  systematic  and  reliable  method  for  the  description  of  the  spin  distribution  within 
large  molecular  systems.  For  radicals  whose  total  electronic  spin  angular  momentum  S 
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equals  or  in  atomic  units'  ^,  interactions  such  as  electron  spin-spin  coupling  do  not 
affect  the  experimentally  measured  spin  distribution,  and  thus  do  not  need  to  be  addressed 
here  and  will  contribute  no  fog  to  the  overview  of  our  final  results.  The  vast  majority  of 
molecular  systems  have  S  values  of  zero;  by  far  the  next  most  prevalent  value  for  S  is 
^.  Seen  in  this  light  then,  our  stressing  spin  one-half  systems  is  not  as  constraining  as 
it  might  seem  at  first,  and  if  our  attempt  at  developing  reliable  computational  methods  is 
successful  then  such  multi-electron  interactions  can  be  addressed  at  our  leisure. 

Let  us  start  with  the  simplest  molecular  ESR  experiment.  Consider  a  molecular 
system  with  a  nonzero  net  electronic  spin.  Associated  with  this  spin  will  be  a  magnetic 
moment.  If  a  static  external  magnetic  field  is  applied,  the  spin  magnetic  moment  of  the 
electron  will  interact  with  it.  If  the  total  electronic  spin  angular  momentum  S  of  the 
molecular  system  is  ^,  and  the  molecular  system  is  spatially  non-degenerate,  then  there 
are  two  magnetic  energy  levels  available  to  the  electron.  These  levels  are  degenerate 
in  the  absence  of  a  second  magnetic  field.  The  degeneracy  is  lifted  by  the  application 
of  the  static  external  magnetic  field  and  the  splitting  is  proportional  to  the  external  field 
strength  (Hq).  The  basic  ESR  experiment  then  consists  of  observing  transitions  between 
these  two  levels,  or  states. 

1.  In  atomic  units  h  =  -e  =  rrie  =  1,  where  h  is  Planck's  constant  divided  by  27r,  e  is  the 
magnitude  of  the  electron  charge  and  is  the  electron  rest  mass.  Unless  otherwise  stated  atomic 
units  will  be  used  through  out  the  rest  of  this  work. 
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Transitions  between  the  energy  levels  can  be  induced  by  the  application  of  a  second, 
appropriately  oriented,  external  oscillating  magnetic  field.  In  our  simple  experiment  the 
oscillating  field  is  applied  at  a  right  angle  to  the  static  external  field  Hq.  Transitions  will 
occur  when  the  resonance  condition  AE  =  hvQ  =  ge^eHo  is  satisfied.  Figure  (1)  below. 


Figure  1:  The  electron  magnetic  resonance  condition. 

Here  h  is  Planck's  constant  (6.63  x  10"^^  erg  sec),  /3e  is  the  Bohr  magneton  (9.27  x 
10'^^  erg/G),  ge  is  a  dimensionless  constant  of  proportionality  equal  to  2.(X)2322  for  a 
free  electron  and  i/q  is  the  frequency  of  the  perpendicularly  applied  oscillating  magnetic 
field.  The  oscillating  magnetic  field  has  the  same  probability  of  inducing  an  upward 
transition  as  a  downward  transition.  In  a  bulk  sample  at  equilibrium  in  the  static  field  the 
lower  level  will  be  more  highly  populated  than  the  upper  one,  thus  there  will  be  more 
upward  than  downward  transitions,  and  a  net  energy  transfer  to  the  bulk  sample.  It  is 
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this  absoqjtion  of  energy  that  is  detected,  or  observed,  experimentally. 

If  this  was  all  there  is  to  an  ESR  experiment  then  little  useful  information  would 
be  obtained  by  its  being  carried  out.  We  would  be  able  to  obtain  the  constant  of 
proportionality  g,  which  would  equal  the  ratio  j^-  If  for  the  moment  we  disregard 
the  possible  dynamic  information  intrinsic  in  the  absorption  line  width  and  shape,  then 
this  would  be  the  only  piece  of  information  we  could  extract.  Fortunately,  this  simple 
picture  is  not  complete.  The  simple  system  we  used  as  our  example  above  is  not  to  be 
the  case  for  the  vast  majority  of  magnetic  molecular  systems.  Almost  certainly  there 
will  be  other  intrinsic  magnetic  fields  within  the  molecular  system.  These  can  be  due  to 
other  unpaired  electrons  or  orbitals  which  possess  a  net  angular  momentum  and  probably 
most  importantly,  to  magnetic  nuclei  within  the  system.  The  interaction  between  the 
unpaired  electron  and  the  various  magnetic  nuclei^  gives  a  wealth  of  information  on  the 
distribution  of  the  unpaired  spin  within  the  system;  this  in  turn  yields  detailed  information 
on  such  varied  subjects  as  the  molecular  geometry,  the  electronic  state  of  the  system, 
the  net  spin  of  the  magnetic  nuclei  can  be  determined,  and  from  this  the  isotope  of  the 
magnetic  center  being  observed.  The  experimentalist  can  also  detect  from  the  electron- 
nuclear  interactions  the  presence  within  the  system  of  nuclear  spin-spin  coupling  as  well 
as  a  wealth  of  other  information  we  need  not  detail  here. 

2.  A  rather  complete  list  of  the  magnetic  nuclei  and  their  various  spin  values  can  be  found  in 
Appendix  (1)  of  Reference  6. 
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From  this  brief  overview  it  is  obvious  that  magnetic  resonance  spectroscopy,  for  this 
work  ESR,  is  a  powerful  tool  for  sciences'  use  in  the  investigation  of  the  invisible  world 
of  atoms  and  molecules.  It  is  just  this  power  to  yield,  sometimes  overwhelmingly  so, 
detailed  and  complex  information  on  the  system  being  observed  that  spawned  the  need 
for,  and  our  desire  to  develop,  theoretical  methods  to  aid  in  explaining,  and  hopefully 
predicting,  the  observed  spectrum. 


The  approximate  molecular  Hamiltonian^  is  a  thirty-five  term  construct,  some  twenty 
of  these  terms  will  influence  the  experientially  observed  ESR  spectrum.  This  dissertation 
shall  primarily  be  concerned  with  only  one  of  these,  the  isotropic  hyperfine  interaction, 
or  the  Fermi  contact  [8],  a,so.  The  origins  of  the  Fermi  contact  term  is  explained  in  terms 
of  a  magnetic  dipole-dipole  interaction  between  the  magnetic  moments  of  an  electron  and 
a  nucleus.  The  vector  potential  [9]  due  to  the  nuclear  magnetic  dipole  moment  is 


The  Isotropic  Hyperfine  Interaction 


(1.1) 


where  we  use  the  standard  mks  units.  The  magnetic  field  due  to  fl^  is 


5(r)  =  V  X  A(r)  =  -  ^- 


(1.2) 


The  interaction  energy  between  jiff  and  /?e  is 


3r  (f  ■  /r^) 


f^N  ■  f^e  _  3(r  •  /7e)  (r  •  /T^) 

|.3  ^5 


(1.3) 


3.  See  Appendix  (F)  of  Reference  7. 
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where  f  represents  a  vector  joining  the  electron  and  nuclear  centers,  /Te  and  fl^  are, 
respectively,  the  electron  and  nuclear  magnetic  moments.  The  Hamiltonian  for  this 
interaction  is  written  as 


Hdip  =  -9e^e9N^N 


S-I^  3(5-r)(/^-r) 


(1.4) 


where  the  collection  of  constants  preceding  the  operators  are,  respectively,  the  free 
electron  "g-factor"  or  Lande  g-factor  ge,  the  Bohr  magneton  ^e.  the  nuclear  g-factor 
and  the  nuclear  magneton  0^-  The  operators  5"  and  /  are  the  total  angular  momentum 
operators  for  the  electronic  system  and  the  nucleus,  respectively. 

When  the  electron  is  found  in  an  orbital  with  an  angular  momentum  of  zero,  i.e. 
/  =  0  or  an  s  state,  the  wavefunction  is  not  zero  at  the  nucleus.  (In  the  Dirac  theory  [4] 
the  atomic  radial  function  diverges,  but  the  divergence  is  so  small  that  equation  (D.7) 
is  approximately  correct  for  the  nuclear  dimension.)  In  any  case  ^,  or  ^,  of  equation 
(1.3)  diverges  at  the  origin.  •  / 

Suppose  we  take  a  small  sphere  with  the  nucleus  as  the  center.  Then  the  contribution 
to  the  isotropic  interaction  from  the  s  electron  outside  the  this  sphere  is  zero,  since  the 
integral  of  p-  for  this  volume  stays  finite.  If  the  sphere  is  small  enough,  the  electron 
density  for  the  s  electron  inside  the  sphere  can  be  regarded  as  a  constant,  and  be  replaced 
by  the  wavefunction  at  the  origin,  i.e.  as  |V'(0)p.  Since  the  s  electron  has  a  magnetic 


moment  /?e,  this  small  sphere  is  homogeneously  magnetized  with  a  magnetization  of 


(1.5) 


The  magnetic  field,  Hs,  at  the  center  of  such  a  homogeneously  magnetized  sphere  is 


M 


(1.6) 


or 


B  =  4 


IT 


(1.7) 


at  the  nucleus.  The  interaction  energy  between  this  field  and  the  nuclear  magnetic  moment 


IS 


(1.8) 


(This  is  the  Fermi-Segr6  term  [10]  and  was  first  given  by  Fermi  [8].)  In  operator  form 
equation  (1.8)  is  written  as 


(1.9) 


where  a,so  equals 


StT  o     Zf  \ 

O-iso  =  -^gePegsPsP  (Tn) 


(1.10) 


Here  />^(r;^)  is  the  unpaired  spin  density  at  the  center  of  nucleus  (N).  As  we  discussed 
earlier,  an  unpaired  electron  in  the  presence  of  a  magnetic  field  has  two  states  available 
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to  it,  either  aligned  with  or  against  the  field.  These  states  usually  are  labeled  spin  up  and 
spin  down  or  as  alpha  (a)  spin  and  beta  spin.  The  unpaired  spin  density  at  a  point 
in  space,  p^{r  —  Vo),  is  then  defined  as  the  difiference  between  the  a  and  ^  spin  densities 
[11]  at  that  point,  or  as  it  applies  to  equation  (1.10)  above 

=  p«(r,)  - /K)  (1.11) 

The  product  ^^^gN^s  is  a  constant.  The  constant  of  proportionality  is  not  altogether 
constant.  Nor  can  it  in  general  be  considered  a  scalar  quantity.  It  is  correctly  expressed 
as  a  second  rank  Cartesian  tensor  [12].  Five  terms  contribute  to  the  tensorial  components 
[7].  Generally  the  largest  contribution  to  the  anisotropy  of  is  through  the  one-electron 
spin-orbit  interaction  [13] 

Hso  =  <;L-S  (1.12) 

where  L  is  the  orbital  angular  momentum  operator  and  {<;)  is  the  spin-orbit  coupling 
constant.  The  one-electron  spin-orbit  interaction  is  proportional  to  the  atomic  nuclear 
charge  raised  to  the  fourth  power,  (Z^),  [13]  and  becomes  important  for  the  heavier 
atomic  systems.  The  molecular  systems  to  which  the  methods  developed  within  this 
dissertation  are  to  be  applied,  are  themselves  constructed  from  the  first,  second  and 
occasionally  the  third  period  elements.  The  one-electron  spin-orbit  interaction  in  these 
atomic  systems  is  almost  always  quite  small.  The  change  in  the  average  value  of  ge 


for  these  systems  is  typically  at  the  third  decimal  place,  an  error  smaller  than  the  error 
intrinsic  to  the  computational  methods  we  will  use  to  determine  the  spin  density.  Thus, 
for  the  methods  developed  within  this  dissertation  the  electronic  g-factor  will  be  treated 
as  a  scalar  constant.  This  leaves  . 

,       _  (1.13) 

as  the  quantity  to  be  determined  computationally. 

Before  describing  the  general  details  of  how  we  go  about  computing  equation  (1.13) 
we  need  to  present  an  overview  of  the  quantum  chemical  approach  we  have  chosen  to 
use.  Our  stated  objective  is  to  derive  a  method  that  can,  to  predictable  limits,  determine 
the  unpaired  spin  distribution  throughout  biologically  important,  and  usually  quite  large, 
molecular  systems,  and  to  test  it  through  the  calculation  of  an  observable  property,  a,so. 
This  constraint  as  well  as  the  desire  to  start  out  on  known  theoretical  ground  leads  one  to, 
as  least  in  my  opinion,  semiempirical  molecular  orbital  (MO)  methods.  Semiempirical 
MO  techniques  such  as  the  various  one-electron  Hiickel  methods  [14,  15]  or  the  two- 
electron  ZDO  [16]  approaches  have  been  applied  to  such  systems  in  the  past  [17,  18]  with 
mixed  results.  Their  simple  form  and  low  computational  cost  make  Hiickel  MO  methods 
logical  first  choices.  Collectively  they  do  an  adequate  job  of  describing  the  gross  overall 
spin  distribution  of  systems  such  as  the  porphyrins  or  the  chlorophylls,  systems  that  shall 
concern  us  later.  On  the  other  hand  the  calculation  of  quantitatively  acceptable  results 
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by  such  methods  has  not  been  demonstrated.  The  assignment  of  the  isotropic  hyperfine 
splittings  (hfs)  through  the  use  of  the  McConnell  relation  [1]  and  Hiickel  spin  densities 
have  proven  to  be  qualitatively  useful.  We  expect  to  do  better  than  this  and  obtain  results 
that  can  be  quantitatively  useful.  The  aforementioned  semiempirical  Hiickel  methods 
do  not  treat  the  electronic  two-body  interactions  of  a  molecular  system  explicidy.  The 
ZDO  semiempirical  methods  we  choose  as  our  theoretical  starting  point  do.  It  would 
seem  from  the  start  something  of  importance  is  missing  from  the  one-electron  Hiickel 
methods.  Also,  primarily  for  the  reason  of  the  ignored  electronic  two-body  interactions 
the  Hiickel  methods  do  not  lend  themselves  to  post  SCF  procedures,  techniques  we  will, 
in  later  chapters,  find  to  be  of  importance  in  obtaining  reasonable  values  for  the  spin 
density,  (1.13). 

Overview  of  Self-Consistant  Field  Theory 

Today  the  ab-initio^  molecular  orbital  theory  of  choice  is  the  Hartree-Fock  Self- 
Consistent  Field  method  [19-25]  (HF-SCF).  The  ZDO  theories  we  have  chosen  to  use  as 
our  starting  point  are  themselves  of  the  Hartree-Fock  self-consistent  field  type.  As  such, 
it  is  necessary  at  this  point  to  outline  Hartree-Fock  theory  before  introducing  the  ZDO 
methods.  It  is  not  our  intent  to  develop  the  Hartree-Fock  equations  in  detail  as  this  lies 
outside  the  scope  of  this  work.  We  do  however  need  to  give  sufficient  detail  here  so 

4.  Ab-initio  theories  tend  to  make  only  those  approximations  which  are  necessary  to  keep  the 
problem  soluble. 
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the  ground  work  will  have  been  laid  for  the  methods  discussed  in  later  chapters.  These 
methods  will  address  the  open-shell  problem  direcdy.  Here  we  will  give  an  overview 
of  closed-shell  Hartree-Fock  theory,  and  leave  the  explanation  of  the  various  open-shell 
approaches  to  each  specific  chapter. 

What  we  would  like  to  be  able  to  solve  is  the  nonrelativistic,  time-independent, 
many-electron  Schrodinger  equation, 

=  (1.14) 

for  the  wave  function  ^,  where  H  is  the  Hamiltonian  for  the  system  of  electrons  and 
nuclei  and  E  is  the  energy  of  the  system.  Explicitly  H  is  written  as 

«=1  A=\        ^  i=\  A=\  i=\  i<3     '  A=1A<B 

Above  Ra  and  are  the  position  vector  and  charge  of  nucleus  A,  and  r.,  the  position 
vector  of  electron  i.  The  summations  run  over  all  N  electrons  and  M  nuclei.  Term 
by  term  equation  (1.15)  is  explained  as  follows,  the  kinetic  energy  operator  for  the  N 
electrons,  the  kinetic  energy  operator  for  the  M  nuclei,  the  Coulombic  nuclear-electronic 
attractions,  the  electron-electron  Coulombic  repulsions  and  lastly  the  Coulombic  nuclear- 
nuclear  repulsions. 

With  two  exceptions  the  wave  function  ^  of  equation  (1.14)  cannot  be  solved  for 
in  closed  form.  These  exceptions  are  the  special  cases  of  1)  the  free  particle  and  2)  a 
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two  interacting  body  system,  or  a  system  that  can  be  reduced  to  two  interacting  bodies. 
The  problem  lies  in  that  the  motion  of  the  electrons  and  nuclei  cannot  be  uncoupled  for 
a  three  or  more  interacting  body  system.  If  we  make  the  following  two  approximations 
we  can  imcouple  the  motions  of  the  various  particles  and  at  least  within  the  validity  of 
the  approximations  solve  for  the  wave  function  These  approximations  are  known  as 
the  Bom-Oppenheimer  approximation  and  the  Independent  Particle  approximation. 

The  Bom-Oppenheimer  Approximation 

If  we  invoke  the  physical  picture  of  the  electrons  moving  among  the  nuclei  so  much 
more  rapidly  than  the  nuclei  themselves  move  that  we  can  consider  the  nuclei  to  be 
fixed  with  respect  to  the  electrons,  then  we  have  also  invoked  the  Bom-Oppenheimer 
approximation  [26].  This  uncouples  the  motions  of  the  nuclei  from  each  other  as  well 
as  from  the  motion  of  the  electrons.  The  effect  of  this  uncoupling  on  the  form  of 
equation  (1.14)  is  as  follows;  the  wave  function  *(ri,i?^)  now  parametrically,  rather 
than  explicitly,  depends  on  the  nuclear  coordinates      and  we  can  write  it  as 

^{r„R^)  =  ^^^{r,)n{R^)  (1.16) 
where  ^(i?^)  is  the  wave  function  describing  the  now  separate  motions  of  the  nuclei  and 


(1.17) 
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is  deemed  the  electronic  wave  function  describing  the  electrons  in  the  field  of  the  "fixed" 
nuclei.  Treating  the  nuclear  positions  as  fixed  parameters  leads  to 

M  ^ 

and 

M    M  ^  ^ 

P  ^  "  Enucrep  =  static  value  (1-19) 

A=1A<B 

Thus  the  Hamiltonian,  equation  (1.15),  simplifies  to 

i=l  ,=1  A=l  .=1  i<j 


(1.20) 


the  electronic  Hamiltonian,  and  equation  (1.14)  is  rewritten  as 

Hel^el  =  E,l^el  (1.21) 
and  called  "the  electronic  SchrOdinger  equation".     .  • 

Even  with  this  considerable  simpUfication  of  the  physical  picture  the  resulting  eigen- 
value problem,  equation  (1.21),  can  only  be  solved  exactly  for  the  special  case  of  one 
electron  moving  in  the  field  of  the  fixed  nuclei.  Ultimately  we  would  like  to  address 
molecular  systems  which  contain  hundreds  of  electrons  but  for  now  we  would  setde  for 

systems  containing  two  to  a  few.  This  leads  to  the  next  great  simplifying  approximation. 
The  Independent  Particle  Approximation 

We  would  like  to  solve  equation  (1.21)  exactly  for  the  many-electron  case  but  the 
electron-electron  interactions  prevent  us  from  being  able  to  separate  the  variables,  (i.e.  the 
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electronic  coordinates)  and  no  closed  form  solution  can  be  obtained.  However,  equation 
(1.21)  can  be  solved  for  the  one-electron  case  which  leads  us  to  a  crude  approximation. 
We  ignore  the  electron-electron  interactions  and  solve  equation  (1.21)  as  N  one-electron 
problems,  where  N  is  the  number  of  electrons  in  the  molecular  system.  The  total 
Hamiltonian  for  the  A^-electron  system  within  this  independent  particle  approximation 
would  become  ; , 

^  =  ■^1  +  ^2  +  ...'^^  vv-  (1.22) 

with  '  . 


(1.23) 


a  sum  of  one-electron  Hamiltonians,  equation  (1.22).  This  form  of  the  Hamiltonian  leads 
to  a  wave  function  of  the  form  - 

*  =  Xi(rOx2(rO---XN(r;.)  (1.24) 

a  product  of  one-electron  fimctions,  orbitals,  and  the  total  energy  of  the  system  will  be 
the  sum  of  the  orbital  energies. 

E  =  +  •  •  (1.25) 

where  the  {x}  in  equation  (1.24)  are  hydrogen-like  wave  functions  or  an  approximation 
to  them. 
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Assuming  that  the  one-electron  Hamiltonian,  equation  (1.22),  and  the  wave  function, 
equation  (1.24),  are  correct  we  can  now  reintroduce  the  electron-electron  interactions  as 
corrections.  This  is  done  in  an  average  way, 

{im  =  J  Jdr,  dr,  x:(i)  X-U)  j-  Xi{i)  xAJ)  (1-26) 

which  represents  the  repulsion  between  two  electronic  charge  distributions,  |xi(OP 
|Xj(i)p.  Adding  these  averaged  electron-electron  interactions  back  into  the  one-electron 
Hamiltonian  yields  the  two-electron  Hamiltonian. 

Hartree  proposed  the  Self-Consistent  Field  (SCF)  method  in  1928  [19].  He  assumed 
the  repulsion  experienced  by  a  particular  electron  is  a  spherically  symmetric  average  of 
the  potential  due  to  the  -  1  other  electrons.  This  leads  to  the  replacement  of  the  exact 
two-electron  potential  of  equation  (1.27)  with  an  effective  one-electron  potential 

^  =  E^/^'  =  EE-^^w  (1.28) 

;  f  - 

where  Jj{i),  the  Coulomb  operator,  is  defined  as 

Jj{i)=  f  dr,^^^^^^^   -         .  .  (1.29) 

When  Jj  operates  on  another  charge  distribution  it  yields  the  classical  interaction  between 
two  charge  distributions.  It  is  now  possible  to  express  the  total  /^-electron  Hamiltonian 
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as  a  sum  of  effective  one-electron  Hamiltonians  of  the  form 


=  hi  +  v;ff 


(1.30) 


leading  to  the  set  of  pseudo-eigenvalue  equations  which  determine  the  one-electron  wave 
functions  {<f>i}. 


Hartree  then  assumed  a  set  of  wave  functions,  orbitals,  and  proceeded  to  calculate  the 
Hartree  operator  for  the  i-th  electron  by  evaluating  the  integrals  of  the  effective  potential. 
The  resulting  N  equations  are  solved  for  new  orbitals  which  are  then  used  to  compute 
the  Hartree  operator  for  the  next  iteration.  This  process  is  considered  finished  when  the 
i  +  1  iteration  is  within  a  given  tolerance  of  the  i-th  iteration  at  which  time  the  potential 
is  termed  self-consistent.  Hartree' s  SCF  method  failed  [27,  28],  this  was  due  primarily 
to  the  form  of  the  wave  function. 

It  is  a  fundamental  principle  of  quantum  mechanics  that  the  measurement  of  any 
physical  property  which  depends  on  the  coordinates  of  a  system  of  identical  particles 
should  not  depend  in  any  way  on  the  arbitrary  labeling  of  the  particles.  For  two  electrons 
described  by  the  wave  function  ^(  x.(l)  X>(2)  )  there  exists  a  probability  density 


(1.31) 


(1.32) 


which  cannot  depend  upon  the  permutation  of  the  electron  labels. 


I*(X.(1)X.(2))|'  =  |*(X.(2)X,(1))P 


(1.33) 
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Which  leads  to  two  possibilities 

x.(l)  X;(2)  )  =      X.(2)  X.(l)  )  (1.34) 

or 

X.(l)  X.(2)  )  =        X.(2)  X.(l)  )  (1-35) 

the  symmetric,  equation  (1.34),  and  antisymmetric,  equation  (1.35),  cases  respectively. 
Experimental  evidence  [29]  indicates  it  is  the  antisymmetric^  case  to  which  the  wave 
function  for  two  electrons  must  conform.  Hartree's  simple  product  form  for  the  wave 
function  is  symmetric  with  respect  to  coordinate  interchange  and  this  is  the  principal 
reason  for  the  failure  of  the  Hartree  method. 

Working  independently  of  one  another  both  Fock  [21]  and  Slater  [22]  applied  the 
variation  principle  to  Hartree's  method.  By  varying  the  radial  functions  of  the  simple 
Hartree  product  to  minimize  the  total  energy  they  rederived  Hartree's  original  intuitive 
method  on  a  more  rigorous  footing,  and  assured  the  total  energy  so  obtained  is  an  upper 
bound  to  the  exact  ground  state  energy. 

When  the  additional  constraints  of  maintaining  orbital  orthonormality  during  the 
variation,  along  with  that  of  an  antisymmetric  wave  function  are  applied  to  Hartree's 
method,  the  Hartree-Fock  equations  are  obtained.    The  constrained  minimization  is 

5.  It  is  a  postulate  of  quantum  mechanics  that  the  wave  function  describing  systems  of  indis- 
tinguishable particles,  which  themselves  possess  half  integer  spin,  must  be  antisymmetric  with 
respect  to  interchange  of  particle  labels. 
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usually  carried  out  using  Lagrange's  method  of  undetermined  multipliers.  While  the 
antisymmetric  wave  fimction  is  constructed  by  the  application  of  A  on  the  simple  Hartree 
product.  The  operator  A  is  the  antisymmetrizer. 

■  P 

The  summation  runs  over  all  A^!  possible  permutations  of  the  particles,  p  is  the  parity  of 
the  permutation  and  is  equal  to  the  number  of  two-fold  permutations  needed  to  bring  the 
orbital  product  back  to  normal  order.  It  was  Slater  who  first  noted  that  equation  (1.36)  was 
the  prescription  for  the  expansion  of  a  determinant  [30].  This  leads  to  wave  functions 
constructed  from  antisymmetrized  products  of  orbitals  being  commonly  referred  to  as 
Slater  determinants. 

The  Hartree-Fock  equations  are  given  as 

N 

ni)M^)  =  Yl^-bM^)  (1-37) 

6=1 

where  /  is  the  Fock  operator 

N  1 


/(I)  Ml) 


Mi)  +  E(W)-^6(i)) 


6=1 


<^a(l)  (1.38) 


The  operator  J  of  equation  (1.38)  is  the  same  Coulomb  operator  as  that  of  Hartree's 
method.  The  operator  K  is  the  exchange  operator  and  gets  its  name  from  the  exchange  of 
electronic  coordinates  that  takes  place  during  its  operation.  The  exchange  interaction  in 
the  Hartree-Fock  equations  arises  due  to  the  antisymmetric  nature  of  the  electronic  wave 
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function,  and  so  did  not  occur  in  the  Hartree  method.  The  Coulomb  operator  was  given 
previously  as  equation  (1.29).  The  exchange  operator  is  defined  as 

^  .      Ki{l)  Ml)  =  y  dr,  m)      <^a(2)  hil)  (1.39) 

Equation  (1.37)  can  be  brought  into  the  form  of  a  standard  eigenvalue  problem  by  a 
unitary  transform  of  the  matrix  of  Lagrangian  multipliers. 

/(I)  <^„(1)  =  £a<^a(l)  ■  '  '  (1.40) 

The  total  electronic  energy  as  obtained  from  the  Hartree-Fock  method  is  given  by 

N  N  N 

I  i  j>i 

where  |*o)  is  the  Slater  determinant  formed  from  the  N  occupied  Hartree-Fock  orbitals. 

Originally  the  Hartree-Fock  method  was  applied  to  atoms  and  the  resulting  differential 
equations  were  solved  numerically.  For  molecular  systems  this  is  impractical  for  all 
but  diatomics.  The  difficultly  arises  in  the  calculation  of  the  Hartree-Fock  orbitals  or 
Molecular  Orbitals.  Roothaan  [23]  and  Hall  [24]  both  suggested  that  the  Hartree-Fock 
orbitals  be  expanded  as  a  linear  combination  of  some  known  set  of  basis  functions 
according  to 

where  is  a  fixed  atomic  orbital  (AO)  basis  function  and  C^i  is  the  expansion  coefficient 
for  the  i/-th  atomic  orbital  in  tiie  i-th  molecular  orbital.    This,  Linear  Combination 
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of  Atomic  Orbitals  -  Molecular  Orbital  (LCAO-MO),  linear  expansion  reduces  the 
numerical  solution  of  the  Hartree-Fock  equations  to  an  algebraic  one.  The  Hartree- 
Fock  problem  becomes  a  system  of  linear  equations  which  are  solved  for  the  expansion 
coefficients  C^i-  The  LCAO-MO  expansion  has  the  added  advantage  of  providing 
molecular  orbitals  with  a  definite  analytic  form,  a  huge  advantage  over  numerical  orbitals 
when  calculating  properties. 

Substituting  equation  (1.42)  into  equation  (1.40)  yields 

V  V 

left  multiplying  (1.43)  by  xjl  and  integrating  over  all  electronic  coordinates  yields 

^t^^^^i  =      SiipC^iei  (1.44) 

V  V 

where 

F^v  =  J  dT.xlil)  fXuil)  (1.45) 

V  =  / dr.xlil)  X.(l)  (1.46) 
If  there  are  atomic  basis  functions  then  there  will  be  molecular  orbitals  and  N 
equations  of  the  form  (1.43).  In  matrix  notation  equation  (1.44)  is  written  as 

FC  =  SCE  (1.47) 

and  is  often  referred  to  as  tiie  Roothaan-Hall  equation.  This  is  the  equation  we  solve  in 
an  iterative  manner  until  self-consistency  is  achieved.  Upon  convergence  the  columns  of 
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C  are  the  coefficients  for  the  AO  expansion  of  equation  (1.42),  one  column  or  MO  per 
original  AO,  and  the  eigenvalues  E  are  the  Hartree-Fock  MO  energies. 

At  this  point  we  need  to  express  the  closed-shell  spin  Restricted  Hartree-Fock  (RHF) 
equations  in  their  fine  detail.  This  will  allow  for  a  better  understanding  of  the  ZDO  based 
approximations  to  them.  Also,  up  to  now  we  have  been  cheating  a  little.  That  is,  we 
have  not  addressed  electron  spin. 

In  1925  Uhlenbeck  and  Goudsmit  [31,  32]  proposed  the  electron  has  an  intrinsic 
angular  momentimi,  or  "spin",  to  explain  the  results  of  the  Stem  and  Gerlach  experiment 
[33,  34].  Therefore  our  description  of  the  electron  using  three  spatial  coordinates  is 
incomplete,  a  forth  coordinate,  the  spin,^  is  needed.  It  was  postulated  that  the  magnimde 
of  this  spin  was  ^  unit  leading  to  eigenvalues  for  Sz  equaling  It  is  customary  to 
designate  this  spin  coordinate  u  and  to  introduce  two  orthogonal  spin  functions  a{uj) 
and  /3{u;)  such  that 

{a\a)  =  W)  =  l  (1.48) 

{a\^)  =  (^|a)  =  0  (1.49) 
5|a)  =  5(5  +  l)|a)  =  f|a)  (1.50) 

m^S{S  +  im  =  i\/3)  (1.51) 

6.  Dirac  placed  the  intrinsic  electron  spin  on  a  solid  physical  footing  with  his  relativistic  treatment 
of  the  electron  in  1928  [4]. 
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Sz\a)  =  l\a)  (1.52) 

5.1/9)  =  (1.53) 

If  we  multiply  the  set  of  molecular  orbitals  {<!>},  or  the  set  of  atomic  orbitals  {x},  by 
these  spin  functions  two  sets  of  "spin"  orbitals  are  obtained.  The  properties  of  the  spin 
orbitals  are  such  that 

tki  =  (1.54) 

^i  =  <l>i\/3)  (1.55) 

(V'.|V'.)  =  (<A.#.)H«)=1  (1.56) 

=  {<i>M{m  =  1  (1.57) 

(V-.l^.)  =           =  0  (1.58) 

or  two  sets  of  spin  orthogonal  orbitals.  How  the  inclusion  of  spin  into  the  wave  function 
affects  its  construction  will  be  explored  in  detail  in  chapters  two  and  three,  equations 
(1.54)  through  (1.58)  above  are  sufficient  detail  for  now.  The  introduction  of  the  spin 
functions  a{uj)  and  /3{u})  into  the  electi-onic  wave  function  has  the  effect  of  enforcing 
tiie  Pauli  Principle.  Two  electi-ons  can  be  described  by,  or  occupy,  one  spatial  orbital 
only  if  they  are  of  different  spin. 

In  Restricted  Hartree-Fock  tiieory  all  occupied  molecular  orbitals  are  doubly  occu- 
pied; there  are  no  half  filled,  or  open-shell,  orbitals  in  RHF  tiieory.  Hence  die  term 
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closed-shell.  If  we  take  our  two  sets  of  spin  orbitals,  and  {ip},  there  will  be  y 
occupied  orbitals  in  each  set,  where  N  is  the  number  of  electrons.  If  we  expand  the  spin 
molecular  orbital  basis  as  was  done  in  equation  (1.42),  this  time  using  the  spin  atomic 
orbital  basis,  again  substituting  the  constructs  into  equation  (1.40),  and  integrate  over  the 
spin  coordinates,  the  restricted  closed-shall  Hartree-Fock  operator  is  obtained. 

N/2 

/(I)  =  hil)  +      [2Ja{l)  -  K^il)]  (1.59) 

a 

where 

Ja{l)=  I  dr,<l>:{2)  —  M2)  (1.60) 

Kail)  =  I  dr,  <j>l{2)  ^  U2)  (1.61) 

The  operator  %  is  a  two-fold  permutation  operator,  operating  on  the  electron  labels.  Its 
effect  on  the  simple  orbital  product  (f>i{l)  (t>j{2)  is 

p,3<^.(l)^j(2)  =  ,^i(2)<^,(l)  (1.62) 

Pi2  exchanges  the  electronic  labels  and  the  permutation  operation  takes  place  prior  to  the 
integration.  The  factor  of  two  in  equation  (1.59)  stems  from  the  fact  that  an  electron  in  one 
molecular  orbital  will  have  two  Coulombic  repulsions  and  one  exchange  interaction  with 
the  pair  of  electrons  occupying  a  second,  different,  molecular  orbital.  There  is  exchange 
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only  between  electrons  of  the  same  spin.  The  energy  of  a  RHF  Slater  determinate  is 


(1.63) 


E  =  2^{a\h\a)  +  Y^[2{ah\ah)  -  {ah\ha)] 

=  2^/i„„  +  ^[2j„i-/i:„6] 

a  0,6 

The  /-th  orbital  energy  is  given  as  ^  - 

ei  =  hii  +  YfiJia-Ki,]  (1.64) 

We  are  now  in  a  position  to  expand  the  closed-shell  Fock  operator  matrix  elements 
into  the  atomic  orbital  basis.  It  is  from  this  vantage  point  that  we  will  be  able  to  best 
understand  the  effect  of  the  various  ZDO  approximations  on  Hartree-Fock  theory.  The 
matrix  elements  of  the  closed-shell  Fock  operator  are  given  by 

N/2 

Wlu)  =  {,i\h\l.)  +  Y^[2{fi\Ja\u)  -  {fi\Ka\l^)]  (1.65) 

a 

where  the  Coulomb  matrix  elements  for  the  a-th  molecular  orbital  are 

Walu)  =  /  dr,  <^:(2)  —  <^a(2)|i/(l))  =  (/iG|i/a)  (1.66) 

and  for  the  exchange  operator 

{fi\Ka\iy)  =  (/x(l)|  [  dr,  ct>l{2)  —  <^a(l)k(2))  =  (/.a|az/)  (1.67) 

where  in  equation  (1.67)  %2  has  already  operated.  The  notation  is  as  follows:  n  and 
V  are  atomic  basis  functions  while  a  is  a  molecular  orbital.  Replacing  Ja  and  Ka  in 
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equation  (1.65)  with  equations  (1.66)  and  (1.67)  respectively  gives 

ffiv  =  hfiv  +  ^[2(//a|z/a)  -  {na\av)]  (1.68) 

a 

an  expression  containing  both  atomic  and  molecular  orbitals.  Expanding  ^q.  in  equation 
(1.68)  above,  in  the  manner  of  equation  (1.42)  yields 

Nl2 

=  V  +  E  E  CAaC:j2(/.<7|i/A)  -  {,.a\Xu)] 


(1.69) 


N/2 

=  V + 2  E  E  ^Aac:j(/.^kA)  - -iHXi.)] 

a  aX 

=  V  +  E  ^A(T[WkA)  -  hfi(T\Xl/)] 

where  the  Pa<t  are  the  first  order  density  matrix  elements. 

N/2 

PXa^^Y^CxaCU  (1.70) 
a 

We  have  now  arrived  at  the  point  where  we  can  discuss  what  the  ZDO  approximation  is 
and  which  of  its  manifestations  we  shall  use. 

i  '  .  i  ■ 

The  INDO  Approximation     "       '  ' 
The  Zero  Differential  Overlap  approximation  is  expressed  as 

XfiXudr  =  0,       /  1/  (1.71) 
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In  words,  for  X/x  7^  Xf  the  differential  in  equation  (1.71)  is  defined  to  equal  zero. 
Thus  the  matrix  of  atomic  overlap  integrals  becomes  the  unit  matrix.  The  effect  of 
this  approximation  on  the  one-  and  two-electron  integrals,  in  terms  of  the  number  of 
such  integrals  and  the  physics  they  contribute  to  the  theory,  is  pronounced.  All  three  and 
four  center  Coulomb  integrals  are  removed,  and  no  exchange  integrals  remain,  or 

^liAi'B  ~  ^/if  (1-72) 

{^^A>^c\l^B(^D)  =  6ft„6xa{fiX\iJiX)  (1.73) 

Throughout  this  introduction  we  have  made  reference  to  the  "various"  ZDO  methods, 
we  are  now  at  a  point  where  we  can  at  least  list  them,  and  outline  the  one  pertinent 
to  this  work. 

Today  there  are  numerous  ZDO  methods  in  use.  They  are  all  in  one  way  or 
another  based  on  of  the  following  approximations,  the  Complete  Neglect  of  Differential 
Overlap  (CNDO)  approximation  [35],  the  Intermediate  Neglect  Differential  Overlap 
(INDO)  approximation  [36],  and  the  Neglect  of  Diatomic  Differential  Overlap  (NDDO) 
approximation  [35].  The  CNDO  method  is  true  to  the  ZDO  approximation,  while  the 
INDO  and  NDDO  methods  relax  somewhat  the  ZDO  constraint  on  the  two-electron 
atomic  integrals.  All  SCF  calculations  performed  within  this  thesis  were  carried  out  at 
the  INDO  level  of  approximation  and  so  it  is  the  INDO  method  we  will  outline  below. 
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Within  the  INDO  approximation  the  Roothaan-Hall  equation,  relation  (1.47),  is 
constructed  and  iteratively  solved  as  was  previously  outlined.  The  difference  between 
INDO  and  ab-initio  Hartree-Fock  lies  in  the  type  of  integral  included,  and  their  functional 
form,  in  the  construction  of  the  Fock  matrix,  with  two  caveats.  The  first  being  that  the 
INDO  method  described  here  explicitly  uses  only  the  valence  atomic  orbitals,  and  the 
second  being  that  in  the  LCAO-MO  expansion  only  one  basis  function  is  used  for  each 
atomic  orbital.  Thus  we  will  be  describing,  and  working  with,  a  valence  orbital  minimum 
basis  set  INDO  method.  Below  the  INDO  equations  analogous  to  equation  (1.69)  are 

presented,  //,  v,  a  and  A  again  label  the  atomic  basis  functions,  {x},  while  A,B,C,  , 

label  the  atomic  center  the  AO  basis  functions  are  associated  with.  The  imposition  of 
the  INDO  approximation  on  equation  (1.69)  leads  to  different  forms  for  specific  cases, 
these  are  listed  below 

1.  Case  fi  =  v;  A  =  B: 

FfiAliA  =  t^/t^/lA  +  [Pec  IaC  -  Vj,c\ 

(1.74) 

Xa  <^A 

2.  Case  n  ^  u;  A  —  B: 

^/ix/ix  =  J]      ^"^^y^  [(/^aA^ Ii^x^^a)  -  \{iiaK  Wa^a)]  (1.75) 

\a  <^a 
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3.  Case  \i  ^  v\  A  ^  B\ 

^t^AVB   —  l^tlAfB  '       2  ^HAVbIaB  (1-76) 

Several  new  terms  were  used  in  expressing  equations  (1.74)  through  (1.76),  these  are 
the  terms  used  throughout  the  literature.  We  shall  explain  them  below  starting  with  the 
one-center  core  integral  (/^^^^  of  equation  (1.74), 


C^Mx/*x  =  (^x|-iv2  +  yj/z^)  (1.77) 


the  one-center  electron-nuclear  attraction, 


V^  =  —  (1.78) 


the  first  order  atomic  density  Pec  is  defined  as, 

Pcc  =  Y.P^cXc  (1.79) 

Ac 

the  two-center  two-electron       is  defined  as, 

Iac  =  (/^A  Ac  1/^.1  Ac)  (1.80) 
and  Vac  tiie  two-center  electron-nuclear  attraction  as 

Vac  =  {iIaWcIiIa)  (1.81) 


29 

The  I3fi^„g  term  is  explained  as  replacing 

^^Ai^B  =  if^A  -        -V,-  V,\i^g)  (1.82) 

A  book  could  be  written,  and  many  have  been  [16,  37-39],  on  the  subject  of  obtaining 
the  various  terms  above.  The  details  of  this  subject  lie  outside  the  scope  of  this  thesis.  It 
should  be  said  the  various  integrals  are  replaced  by  functions  of  experimentally  obtained 
parameters  such  as  ionization  potentials  and  electron  affinities,  with  the  exception  of  in 
equation  (1.82).  The  /3  term  is,  in  many  implementations,  a  totally  empirical  parameter 
and  is  usually  obtained  by  varying  it,  while  calculating  a  specific  property  or  reproducing 
results  of  a  higher  level  of  theory,  until  acceptable  results  are  obtained.  Such  is  the 
case  in  this  work,  we  shall  obtain  the  SCF  orbitals  via  two  INDO  parameterizations,  one 
having  been  parameterized  to  reproduce  molecular  structures  and  the  second  to  reproduce 
the  low  energy  ultraviolet  and  visible  absorption  spectra  of  molecular  systems.  Both 
parameterization  follow  the  philosophy  of  Pople's  original  lNDO/1  [36]  method,  at  least 
with  respect  to  the  functional  form  of  the  U  and  7  integrals  and  in  the  sense  of  varying 
/3  to  reproduce  a  benchmark  value.  ' 

We  are  now  in  a  position  to  detail  how  one  actually  evaluates  the  spin  density,  (1.13). 
At  this  point  we  will  limit  the  discussion  to  a  single  Slater  determinant.  The  point  was 
made  earlier  that  the  unpaired  spin  density  at  any  point  in  space  (r)  was  defined  to  be 
the  alpha  density  at  point  (r)  minus  the  beta  density  at  point  (r).  In  terms  of  the  spin 
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molecular  orbitals,  equations  (1.54)  and  (1.55),  the  alpha  and  beta  densities  evaluated 
at  (r)  are  expressed  as 

N" 

P"ir)  =       IC(r)P  (1.83) 

a 

and 

TV" 

/(r)  =  E  \'^a{r)f  (1.84) 

a 

Then  ^^(r)  =  ^"(r)  —  p^{r)  and  if  we  expand  equations  (1.83)  and  (1.84)  in  the  manner 
of  equation  (1.42)  we  obtain 

P^i^)  =  T.T.P>'^(^)x:{r)  (1.85) 
/i  1/ 

and 

^'w^EE^if-x.wxrw  (1.86) 

where  P*^  and  P^v  are  again  elements  of  the  first  order  density  matrices  of  the  alpha 
and  beta  spaces  respectively,  or 

N" 

p^,.  =  Y.^uc^.ay  (1.87) 

a 

-    .       PS.  =  J2^UCSa)\,  (1.88) 
a 

Again  the  Cft„  are  the  expansion  coefficients  obtained  from  the  solution  of  the  Roothaan- 
Hall  equation,  relation  (1.47).  The  expression  for  p^{r)  can  now  be  written  as 

P'i^)  =  E  E         -  PS.)x.{r)  Xlir)  (1.89) 


31 

where  all  quantities  in  equation  (1.89)  above  have  been  previously  defined.  We  should 
note  that  for  the  closed- shell  RHF  theory  we  previously  outlined,  the  term  ^P*^  -  P^^^  is 
everywhere  identically  zero,  a  consequence  of  each  spatial  orbital  being  doubly  occupied. 

We  can  now  proceed  with  the  thrust  of  this  work,  that  is,  to  develop  a  method  capable 
of  quantitatively  describing  the  unpaired  spin  distribution  throughout  a  molecular  system 
as  measured  through  the  calculation  of  the  isotropic  hyperfine  interaction,  aiso,  and  to  be 
able  to  apply  this  method  to  systems  that  are  of  interest  to  chemists  today.  In  particular, 
we  will  eventually  be  interested  in  very  large  systems,  and  so  our  aim  will  be  more  toward 
applying  these  methods  to  semiempirical  quantum  chemistry.  However,  the  theory  we 
develop  should,  in  general,  be  appropriate  for  any  wave  function. 


CHAPTER  2 
THE  UHF  APPROACH 

Introduction 

In  Chapter  1  we  developed  RHF  theory  in  sufficient  detail  as  to  allow  for  the  formal 
definition  of  p^{r).  Unfortunately,  RHF  theory,  by  construction,  cannot  address  molecular 
systems  which  possess  unpaired  electrons.  The  simplest  method  for  addressing  such 
systems  at  the  SCF  level  of  theory  is  the  Unrestricted  Hartree-Fock  (UHF)  method  which 
is  also  known  as  the  Different  Orbitals  for  Different  Spins  (DODS)  [40-45]  approach. 
Historically  this  method  dates  to  the  late  1920s  where  both  Hylleraas  [46]  and  Eckart 
[47]  applied  it  to  the  He  atom,  and  the  method  as  it  is  formulated  today  dates  to  the  early 
1950s.  Within  UHF  theory  electrons  are  still  assigned  to  orbitals  in  Slater  determinants 
[36],  but  the  spatial  orbitals  associated  with  a-spin  electrons  are  allowed  to  differ  from 
those  describing  the  /3-spin  electrons.  This  additional  freedom  has  several  consequences; 
due  to  the  increased  number  of  variational  parameters  the  UHF  open-shell  energy  will 
be  <  the  ROHF^  energy,  as  well  as  allowing  for  a  more  accurate  treatment  of  electron 
correlation  at  the  SCF  level  [48-50].  Perhaps  the  most  important  ramification,  at  least  for 
this  work,  is  that  a  single  Slater  determinant  constructed  from  UHF  spin  molecular  orbitals 

1.  Both  the  Restricted  Open-Shell  Hartree-Fock  (ROHF)  method  and  the  concept  of  electron 
correlation  will  be  treated  in  Chapter  3. 
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will  not  in  general  be  an  eigenfunction  of  [51,  52],  though  it  will  be  an  eigenfunction 
of  Sz-  As  we  shall  see  later  in  this  chapter  this  so  called  "spin  contamination"  will 
prove  to  be  a  problem  in  the  calculation  of  aiso  in  the  larger  radical,  molecular  or  ionic, 
systems.  We  now  need  to  develop  UHF  theory  in  sufficient  detail  to  understand  how 
Puhf{^)  determined. 

The  UHF  Equations  «• 
In  Chapter  1  we  found  the  general  Hartree-Fock  eigenvalue  equation  to  be 

f{l)rPa{l)  ^  SaM^)      .  (2.1) 

where  the  {ipi}  are  the  restricted  spin  molecular  orbitals.  We  shaU  use  equation  (2.1) 
as  our  starting  point  in  developing  UHF  theory.  In  a  manner  analogous  to  that  used  in 
defining  equations  (1.54)  and  (1.55),  the  restricted  spin-orbital  definition,  we  define  two 
sets  of  unrestricted  spin-orbitals. 

^i  =  (f>i\a)  (2.2) 

i^^  =  <l>M  (2.3) 

where  the  restriction  that  (f>^  =  <^f  imposed  in  RHF  or  ROHF  theory  has  been  dropped, 
thus  giving  the  method  its  name.  Unrestricted  Hartree-Fock.  If  we  now  replace  the  {ipi} 
in  equation  (2.1)  with  the  forms  defined  above,  left  multiply  by  (a|  and  (/?|,  respectively. 
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and  integrate  the  resulting  equations  over  the  limits  of  the  spin  variable  u  we  obtain 

r(l)         =  el  4>^{\)  (2.4) 


and 


r(i)  <^f(i)  =  £f  <j>i{i)    ^  (2.5) 


where 


-    r  (1)  =  h{i)  +    [-^a  (1)  -      +  E  -^'(1)  (2.6) 
r  (1)  =  Ml)  +  E  (1)  -  ^6(1)]  +  E  (1)  (2.7) 

6  a 

the  a-spin  and  /?-spin  Fock  operators  respectively.  The  summation  upper  limits  p  and 
q  are  respectively  the  number  of  a-spin  and  ^-spin  electrons,  and  by  convention  p  >  q. 
Since  the  kinetic  energy  of  the  electron  and  coulombic  nuclear  attraction  are  independent 
of  spin,  h{l)  is  identical  to  that  defined  in  equation  (1.23).  The  unrestricted  coulomb  and 
exchange  operators  are  defined  below,  the  notation  is  the  same  as  that  used  in  defining 
equations  (1.60)  and  (1.61). 

J^(l)  =  / dr,  ra*{2)  —  ra{2)  (2.8) 

J?(l)  =  /  ^r,      (2)  ^  ^,^(2)  (2.9) 

K^il)  =  /  dr,  <^r(2)  —  <^a(2)  (2.10) 
J  '"12 

^,^(1)=  /  dr,r:{2)^<f>t{2)  (2.11) 

J  '12 


35 


The  energy  of  a  UHF  Slater  determinant  is  given  by 


f  1  r 


p  1  1  p 


a  b  a,b 


(2.12) 


a,b  a.  b 


where  equation  (2.12)  is  analogous  to  that  given  for  the  RHF  energy  of  a  Slater 
determinant,  equation  (1.63).  The  J^j  term  is  defined  as 


the  coulombic  repulsion  between  the  a-spin  and  ^-spin  electrons  whose  spatial  distribu- 
tions are  described  by  \<f>l\^  and  |<?i>jp  respectively,  aU  other  terms  in  equation  (2.12)  are 
analogous  to  their  RHF  counterparts. 

We  now  expand  our  two  sets  of  spin  molecular  orbitals  as  a  linear  combination  of 
some  known  set  of  basis  functions  (i.e.  The  LCAO-MO  approximation.)  consisting  of 
m  linearly  independent  functions  {x}. 


drjr,  <^r(l)C(2)  —  KmiC^)  (2.13) 

^12 


m 


(2.14) 


m 


(2.15) 


If  we  now  substitute  equation  (2. 14)  into  equation  (2.4)  we  obtain 


(2.16) 
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now  multiplying  equation  (2.16)  by  xjl(l)  and  integrating  over  yields 


Efa  /-la   


(2.17) 


V 


V 


where  all  quantities  in  equation  (2.17)  were  defined  for  equation  (1.44),  and 


(2.18) 


In  matrix  notation  equation  (2.17)  is  written  as 


(2.19) 


Similar  results  are  obtained  when  one  starts  by  substituting  equation  (2.15)  into  equation 
(2.5)  and  carrying  out  the  subsequent  integration,  or 


These  two  equations,  (2.19)  and  (2.20),  are  the  unrestricted  generalizations  of  the 
Roothaan-Hall  equation,  (1.47).  They  were  first  given  by  Pople  and  Nesbet  [40],  and  are 
often  referred  to  as  the  Pople-Nesbet  equations. 

Once  the  matrices  C°  and  have  been  determined  for  a  molecular  system,  /9^(r)  is 
determined  from  the  following.  If  A  is  a  rectangular  matrix  of  order  m  x  p  constructed 
from  the  first  p  singly  occupied  orbitals,  columns,  of  C",  and  B  is  a  rectangular  matrix 
of  order  m  x  9  constructed  from  the  first  q  singly  occupied  orbitals,  columns,  of  C,  then 


(2.20) 


(2.21) 
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P"  =  BB^  (2.22) 

the  a-density  and  ^-density  matrices  respectively,  which  are  used  to  construct  the  total 
density  matrix 

=  P"^  +  P^  (2.23) 

which  follows  directly  from  equations  (1.87  -  1.89)  and  uiso  is  obtained  from  equation 
(1.10). 

Spin  Contamination 

While  the  UHF  method  has  the  important  advantage  of  directly  producing  and 
P^^  at  the  SCF  level,  it  also  has  several  drawbacks  as  far  as  this  work  is  concerned.  It 
suffers  from  the  same  problem  any  single  determinant  method  has  in  describing  an  open- 
shell  system  which  is  not  of  the  highest  multiplicity,  (e.g.  A  singlet  system  consisting 
of  two  open-shell  orbitals  describing  two  electrons.)  This  is  not  an  insurmountable 
problem  as  the  molecular  systems  in  which  we  are  interested  rarely  possess  sufficiently 
high  symmetry  for  there  to  be  non-accidental  orbital  degeneracies;  of  course,  accidental 
degeneracies  or  near  degeneracies  will  be  a  problem  when  they  occur.  The  crucial 
problem  lies  in  how  "^uhf,  not  in  general  being  an  eigenfimction  of  S\  affects  PuHri^)- 
As  we  shall  see  later  the  affect  on  Puhf{^)  can  be  sufficiently  pronounced  so  as  to  render 
the  calculated  a,so  both  quantitatively  and  qualitatively  useless.  For  now  we  need  to 


■( 
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present  an  overview  on  the  origins  of  spin  contamination  so  the  sections  dealing  with 
minimizing  its  effects  upon,  or  removing  it  from,  Pujjp{r)  will  be  more  lucid. 
A  single  determinant  UHF  wavefunction  is  of  the  form 

=  ^|<^^(l)a(l)  cf>;{p)a{p)  <^f (p  +  l)^{p  +  1)  <(n)/3(n)|  (2.24) 

where  in  the  basis  of  UHF  canonical  spin  molecular  orbitals 

imj)  =  im',) = (2.25) 

and  in  general 

mr,)^o  (2.26) 

The  functions  {(j)"}  and  {<f>^}  can  be  transformed  among  themselves  with  unitary  trans- 
forms 

Vi=J2rsUsi  (2.27) 

Vi^Y^</>lVsi             />  (2.28) 

and  still  leave  (2.24)  unaltered.^  , 
Now  by  applying 

S'  =  S-S+  +  Sl  +  S,  (2.29) 

2.  See  Appendix  A. 
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to  (2.24),  it  is  readily  shown  that 

p  9 

{^vhAS'I^uhf)  ^Sl  +  S,+q-Y,Yl (2.30) 

«  j 

Equation  (2.30)  is  general  and  holds  for  any  single  determinant  wavefunction.  (e.g.  RHF, 
ROHF  as  well  as  UHF.)  It  is  immediately  clear  that  for  RHF  and  ROHF  wavefunctions 
the  last  two  terms  on  the  right-hand-side  (rhs)  of  equation  (2.30)  cancel,  (i.e.  The  overlap 
matrix  in  the  basis  of  restricted  Hartree-Fock  molecular  orbitals  is  the  unit  matrix.)  and 

RHF /ROHF  —  5^(5^4-1)  (2.31) 

However  for  the  UHF  wavefunction 

{S')aHF  >  S,iS,  +  1)  (2.32) 

where  the  equality  in  equation  (2.32)  occurs  when  the  UHF  wavefunction  converges  to 
an  exact  eigenfunction  of  S^,  the  RHF/ROHF  solution,  else 

UHF  EXACT  (2.33) 

Other  Interesting  Orbitals 

The  unrestricted  molecular  orbitals  are,  of  course,  not  the  only  interesting  type  of 
orbital.  We  shall  find  two  other  types  of  orbitals  to  be  of  use  in  understanding  spin 
contamination  in  a  UHF  wavefunction  and  in  removing  it.  These  are  known  as  the 
"natural"  orbitals  [43],  which  are  those  orbitals  that  diagonalize  the  reduced  first-order 
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density  matrix;  and  the  "corresponding"  orbitals  of  Amos  and  Hall  [43]  which  are  also 
known  as  the  "paired"  orbitals  of  Lowdin  [53].  These  are  the  orbitals  which  diagonalize 
the  overlap  matrix  between  the  a-spatial  orbitals  and  the  ^-spatial  orbitals. 
For  the  wavef unction  (2.24)  the  first-order  density  matrix  is 

7(1|1')  =  j2  ^r(l')  +  E  '^'^(l)  ^'{l')  r  (l)^(l')  (2.34) 

r  a 

This  is  already  in  diagonal  form  and,  since  the  canonical  spin-orbitals  are  already 
orthonormal,  they  are  the  natural  spin-orbitals.  The  natural  spin-orbitals  are  not  unique 
since  the  unitary  transforms  (2.27)  and  (2.28)  leave  the  form  of  (2.24),  and  (2.35), 
unchanged.^ 

The  reduced  first-order  density  matrix  is  obtained  by  integrating  (2.34)  over  the  full 
range  of  the  spin  functions  a{u)  and  ^{uj),  this  yields 

=^r(l)<^^(l')  +  E'^r(l)«^f(l')  (2.35) 

r  s 

which  would  be  in  diagonal  form  except  that  {(f)"}  and  {4>^}  are  not  mutually  orthogonal. 
The  unitary  transforms  (2.27)  and  (2.28)  can  be  used  to  define  new  orbitals  having  the 
property  that  their  spatial  overlap  matrix  is  diagonal,'^  or 

Trs  =  J  drr^frjl  [  (2.36) 

3.  See  Appendix  A. 

4.  When  the  unitary  matrices  U  and  V,  of  equations  (2.27)  and  (2.28),  are  determined  as  in 
Appendix  C. 
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Tts  =  TrSrs  (2.37) 

these  are  the  "corresponding"  orbitals  of  Amos  and  Hall. 
In  terms  of  the  corresponding  orbitals  (2.35)  becomes 

Kill')  =  Era)  <(!') +  1:^(1)^1(1')  (2.38) 

r  s 

where  the  {7/°}  are  spatially  orthogonal  to  each  other  and  to  the  {tj^}  except  when 
r  —  s.  If  these  sets  of  paired  functions  are  replaced  with  the  Schmidt  orthogonalized 
[54]  functions 

K  =  {r,;  +  r]'^){2  +  2Tr)-'^  (2.39) 

f^r  =  in;  -         -  2Tr)-^  (2.40) 
then  all  the  functions  are  orthonormal  and  (2.38)  becomes 

/'(i|i')  =  E(i+3^r)A:(i)A,(i')+E(i-^^)/'**(i)/^^(i')+  E  ^r(i)<^t(i') 

(2.41) 

the  functions  A^,  ^ls  and  of  relation  (2.41),  are  the  natural  orbitals.  The  corresponding 
orbitals  are  then  the  natural  spin-orbitals  most  easily  related  to  the  natural  orbitals,  and 
the  overlap  integrals  of  the  "paired"  orbitals  determine  the  natural  orbital  occupation 
numbers. 

A  second  approach  can  be  taken  in  determining  the  UHF  natural  orbitals,  one  that 
will  give  further  insight.  First  we  express  the  reduced  first-order  density  matrix  in  terms 
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of  the  expansion  basis  {x} 

m 

p(l|l')  =Y.(Prs+Pr's)x:{l)Xs{l')  (2.42) 
r,3 

where  P^^s  and  P^^,  are  matrix  elements  of  relations  (2.21)  and  (2.22)  respectively.  To 
see  the  relationship  between  definitions  (2.41)  and  (2.42)  we  first  define  the  matrix 

C  =  [ABU„  (2.43) 

where  [AB]  is  the  concatenated  matrix  formed  by  matrix  A  of  relation  (2.21)  and  matrix 
B  of  relation  (2.22),  mdn  =  p  +  q.  Then  from  equations  (2.21)  and  (2.22)  it  follows  that 

CC^  =  (P'"+P^)„xm         '  (2.44) 
and  from  Appendix  B  this  has  the  same  non-zero  eigenvalues  as  C+C.  Now 

C^C=(h  (2.45) 


where  I,  and  I,  are  the  unit  matrices  of  order  p  x  p  and  q  x  q  respectively,  and  S 
is  the  overlap  matrix  between  the  UHF  a-spatial  orbitals  and  ^-spatial  orbitals.  When 
the  unitary  matrices  U  and  V  are  determined  as  in  Appendix  C,  relation  (2.45)  can  be 
brought  into  the  form 

v^)(s^  0(o  v)=e.  I) 

where  T  is  the  diagonal  matrix  of  corresponding  orbitals  overlaps.  There  are  then  2q 
eigenvalues  of  the  form 

1  ±  Tr  (2.47) 
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and  {p  -  q)  eigenvalues  that  are  exactly  unity,  as  in  (2.41).  If  the  number  of  expansion 
basis  functions  {x}  is  less  than  the  number  of  electrons  in  the  molecular  system,  (i.e. 
m  <  n),  then  there  will  be  m  -  n  overlap  integrals,  equation  (2.36),  that  equal  unity. 
This  corresponds  to  a  set  of  orbitals  whose  occupation  numbers  equals  2.  There  will  be 
n  -  m  =  D  of  these  doubly  occupied  orbitals,  and  they  have  no  pairing  properties.  If 
m  >  n,  there  will  be  a  set  of  orbitals  which  have  occupation  numbers  of  zero.  There 
will  hem-n  =  E  such  orbitals,  and  they  also  have  no  pairing  properties.  When  m  =  n, 
then  all  overlap  integrals  between  the  paired  orbitals  will  be  less  than  unity  and  greater 
than  zero.^  i.e. 


n  >  m,  D 


n  —  m 


(2.48) 


E 


0 


n  <  m 


D 


0 


(2.49) 


E 


m  —  n 


n  =  m 


D  =  E  =  0 


(2.50) 


These  relationships  can  be  stated  clearly  as  the  eigenfunction  equation 


(2.51) 


5.  See  Appendix  C. 
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where  'p  is  defined  as 


P  =  Y,dkOk{\)ei{l')  (2.52) 

*;=1 


then 


Ok=rk  =  ril  dk  =  2^1  +  Tk,     k  =  l,2,-  -D  (2.53) 

ek  =  \k,  dk  =  l  +  Tk,      k  =  D  +  l,---q  (2.54) 

0k  =  tik,  dk  =  l-Tk,      k  =  D  +  l,---q  (2.55) 

ek  =  r]l  =  (f>l,   dk  =  l,          k  =  q  +  l,---p  (2.56) 

There  are  no  relations  for  the  unoccupied  orbitals  {E}. 
Spin  Contamination  and  the  UHF  PE  Surface 

For  RHF  and  ROHF  wavefunctions  the  natural  orbitals  are  the  canonical  molecular 
orbitals  and  their  occupation  numbers  are  either  2,  1  or  0.  It  then  follows  that  for  an 
exact  spin  eigenfunction 

occ 

^aJ  =  49  +  (p-9)  (2.57) 

i 

where  the  o-,-  are  the  natural  orbital  occupation  numbers  [55].  For  a  single  determinant 
no  c7,  other  than  the  RHF  and  ROHF  values  can  satisfy  (2.57)  while  conserving  the  total 
number  of  electrons.  For  a  UHF  wavef unction  the  natural  orbital  occupation  numbers 
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are  functions  of  tiie  corresponding  orbital  overlaps,  or, 

(7i  =  l  +  Ti,     i  =  l,2,--  -  q  j"  ... 

(7.-  =  1-r,',    i'  =  l,2,---  q  (2.58) 

(7^  =  1,  u^2q  +  l,2q  +  2,---  p 

where  i  and  i'  are  paired.  Thus  for  a  single  determinant  UHF  wavefunction  to  be 
an  exact  spin  eigenfunction  tiie  T,  must  equal  1,  i.e.  tiie  corresponding  orbital  overlap 
matrix  must  be  the  same  as  that  of  an  RHF  or  an  ROHF  wavefunction.  This  implies  tiiat 
the  UHF  wavefunction  must  be  related  to  tiie  RHF  or  tiie  ROHF  wavefunction  tiirough 
suitable  unitary  transformations. 

To  this  point  we  have  expressed  the  Fock  operator  as  a  function  of  the  Fock  molecular 
orbitals.  We  could  have  just  as  easily  expressed  it  as  a  function  of  the  total  density  P"^,  let 

/(pT)  =  h  +  J(P^)  -  lK(P'')  (2.59) 
^(P")  =  +  P'^)  /      ^*(2)  ^  u{2)  (2.60) 

K{Fn  =  J2(P^,  +  P'.)  /rfr,/(2)^K2)  (2.61) 

J  '12 

HI/ 

Then  in  the  basis  of  the  corresponding  orbitals  we  can  write  the  UHF  equations  as 

(^r l/(P^)  -  iK{F')\r,-)  =  0  ,  i  =  1, 2,  •  •  •  p,  1/  =  p  +  l,p  +  2,  •  •  •  m  (2.62) 
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{vf\f{Pl  +  lK{F'M)=0,  i  =  l,2,---  9,  u  =  q  +  l,q  +  2,---  m  (2.63) 

where  =  -  P^.  For  an  ROHF  like  solution  the  {r]"}  and  {rj^}  must  be  the  same, 
or,  (2.63)  can  be  written  as 

{Vi\f{P^)  +  ^K{F')\vt)  =  0,i  =  l,2,---  p,  u  =  p+l,p  +  2,---  m  (2.64) 

Comparing  equations  (2.62)  and  (2.64),  it  is  obvious  that  they  both  cannot  be  simulta- 
neously correct  without  P^  =  0,  identically.  This  corresponds  to  an  RHP  solution,  but 
from  relation  (2.57)  an  open-shell  spin-exact  single  determinant  wavefunction  must  be 
related  to  an  ROHF  wavefunction.  Thus  no  single  determinant  spin-exact  wavefunction 
is  stationary  on  the  UHF  potential  energy  surface.  This  result  will  have  some  bearing 
on  the  next  two  sections. 

Spin  Projected  UHF 

As  noted  above  the  single  determinant  wavefunction  (2.24)  is  an  eigenfunction  of  Sz 
with  {Sz)  =  ^{p-q).  However  it  suffers  from  spin  contamination,  described  previously, 
and  for  the  aforementioned  reasons  is  not  an  eigenfunction  of  S'^.  Since  we  use  the 
PUHF  method  to  "improve"  our  UHF  wavefunction,  we  will  describe  below  this  method 
(spin  projection)  for  removing  the  unwanted  spin  contamination  intrinsic  to  the  open-shell 
UHF  wavefunction.  As  we  shall  see  when  computed  properties  are  ultimately  the  goal, 
it  is  far  easier  to  work  with  the  densities  than  it  is  to  work  with  the  wavefunction.  The 
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method  we  detail  below  is  then  centered  on  how  the  densities  are  affected  as  the  spin 
contamination  is  diminished  or  removed.  It  was  not  our  intention  to  rederive  the  algebra 
for  projecting  out  the  5  =  component  of  the  spin  density  matrix.  We  instead  give  an 
overview  and  list  the  pertinent  equations  and  citations. 

When  expressed  in  terms  of  the  corresponding  orbitals  the  UHF  wavefunction  can 
be  written  in  a  form  similar  to  that  of  a  Configuration  Interaction  wavefunction  [56,  57], 
e.g.  equation  (3.15),  or 

where  the  Cs  are  functions  of  the  corresponding  orbital  overlap  integrals.  The  ^  are  linear 
combinations  of  Slater  determinants  constructed  from  the  corresponding  orbitals,  and  are 
eigenfunctions  of  S\  The  excited  terms,  •  •  •,  will  contain  the  5^  =  |(p  -  q) 

componentof  the  pure  spin  states  1^1  =  ^{p  —  q),       — 9  +  |(p-?  +  e/),  where  e/ 

is  the  excitation  level.  It  is  these  higher  multiplicity  states,  l^l  >  ^{p-  q),  that  will  be 
removed  upon  projection.  Relations  for  the  Cg  were  given  by  Sasaki  and  Ohno  [58,  59] 
and  are  reproduced  here  since  they  were  used  in  analyzing  UHF  wavefunctions  in  this 
work.  In  equation  (2.66)  below,  ujs  is  the  "weight"  of  the  components  with  multiplicity 
(2s  +  1),  or,  the  fraction  of  the  wavefunction  composed  from  the  {2s  +  1)  components, 

--S 

Us  =  CI  =  (26  +  l)^^i^  ^(-l)^i?,B3-^.+,-  (2.66) 
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where  s  and  are  the  quantum  numbers  associated  with  and  Sz.  While  s  is  a 
summation  index,       =  ^(p  —  9). 

_[(.-m.  +  i)!p 

i!(2.  +  j  +  l)!  ^  ^ 

(2) 

1=1 

and  the  T,  are  those  of  relation  (2.36).  If  the  wavefunction  is  normalized,  then 

X;^s  =  l  (2.69) 

To  obtain  an  eigenfunction  of  a  projector  [41,  60,  61],  Os,  is  applied  to  (2.24) 
or  (2.65),  where 

ds  =  l[At  (2.70) 


and 


01  =  O 


s 


(2.71) 

S'  Os  =  S{S  +  1)  ds 
is  the  normalized  annihilation  operator  [41,  62] 

AS-      ^^-^(^-1)  .2  72) 

^'  -  ,(3  -  1)  -      -  1)  (2.72) 

and  t,  of  equations  (2.70)  and  (2.72),  ranges  over  all  values,  except  t  ^  s,  consistent 
with  the  number  of  electrons  in  the  system.  The  projection  operator  (2.70)  projects  out 
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the  {2s  +  1)  component  of  (2.65).  The  projected  state  wavefunction  can  be  written  as  a 
linear  combination  of  determinants  [41]. 

O.^^H.  =  ^Ck{s,m,,^)  T^^  (2.73) 
k=o 

The  Ck{s,mz,^)  [58,  59]  are  determined  as 

r(  n,  ,  ,,(f-m,-fc)!(3  +  m,)! 
C7,(.,m„  ^)  =  (2.  +  1)  

(2.74) 

Z.^     ^  /!(s_^^  +  /_A;)!(f-5-/)!(26  +  /+l)! 
The  T^'  are  a  sum  of  Slater  determinants  constructed  from  the  corresponding  orbitals 

[63,  64],  note  that  (2.75)  is  not  normalized. 

Tr  =  I'/r  •  •  •  ri;,  v'^r--  n:,M'-'/3Vl3'-']  (2.75) 

The  [a^~*^*^|  is  the  sum  of  all  spin-function  products,  of  the  p  —  a(u;)  and  k  j3{uj) 
spin-functions,  e.g. 

[a'-' I  =  aa/?  -1-  a^a  +  jSaa  (2.76) 
and  |a*^'~*^]  is  similarly  defined.  Then  [aP~*^*|a''^'~^]  is  a  sum  of 

terms  each  containing  p  a{uj)  factors  and  q  /3{u))  factors  [41,  65]. 

Harriman  [66]  showed  that  the  above  described  projection  leaves  the  eigenfunctions 
of  the  spin  density  matrix,  as  well  as  those  of  the  first-order  reduced  density  matrix,  (2.35), 
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unchanged.  Only  the  eigenvalues  change  upon  projection  and  both  sets,  projected  and 
unprojected,  can  be  expressed  as  functions  of  the  corresponding  orbital  overlap  integrals, 
equation  (2.36),  and  the  {Cjt},  equation  (2.74). 

In  terms  of  the  natural  orbitals  of  the  reduced  first-order  density  matrix,  (or  the  charge 
density  matrix)  relations  (2.39-2.41),  the  natural  orbitals  of  the  spin  density  matrix  are 
expressed  as  -  \ 

<^f  -  ^  (A.  + /^.)  ,  (2.78) 

•^f' =      (A.  -  (2.79) 

<l>t  =  K  (2.80) 

In  eigenfunction  form  the  relationships  are  expressed  as 

roi  =  ikei      '  (2.81) 

where  ^  is  defined  as 

r  =  Y.ikei{i')ei\i)  (2.82) 


and 


n  =  <l>l     6  =  (l-r|)^      k  =  D^\,--q  (2.83) 
ei^<t>l',     ik  =  -{l-Tl)^,     k^D  +  l,---q  (2.84) 
ei  =  <l>%,     6  =  1,     k  =  q  +  l,---p  (2.85) 
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As  was  the  case  for  the  charge  density  natural  orbitals,  there  are  no  relations  for  the 
unoccupied  orbitals  {E}.  The  Tk  of  the  doubly  occupied  orbitals  {D}  all  equal  unity 
and  thus  have  occupation  numbers  of  zero,  see  equations  (2.83)  and  (2.84)  above,  and 
so  do  not  contribute  to  the  spin  density  matrix.  The  validity  of  relations  (2.78-2.80) 
can  be  checked  by  expanding  them  in  terms  of  the  corresponding  orbitals  and  applying 
—      —  p~,  where 

^  =  E^."(l')'7r(l)  (2.86) 

r 

P^  =  E'?f(l')'?r(l)  (2.87) 

s 

and  check  against  relations  (2.83-2.85). 

If  we  now  follow  Harriman's  [66]  line  of  reasoning  we  can  obtain  the  above  results 
from  a  different  perspective.  We  also  will  arrive  at  a  point  where  Harriman's  observations 
on  the  effects  of  projection  are  apparent. 

We  now  want  to  construct  the  operator  matrix  ^  in  the  basis  of  the  charge  density 
natural  orbitals,  and  then  obtain  its  eigenvalues  and  eigenfunctions.  First  we  should  note 
that  can  couple  differing  orbitals  only  if  the  differing  orbitals  are  one  of  the  pairs  of 
"paired"  orbitals,  relations  (2.39)  and  (2.40).  Then  the  largest  non-zero  sub-blocks  will 
at  most  be  of  order  2x2,  and  lie  along  the  principle  diagonal.  This  reduces  the  problem 
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to  diagonalizing  a  series  of  2  x  2  matrices  of  the  form 

f{K\r\^,)  .2  88) 

to  obtain  the  eigenvalues,  and  since  (2.88)  is  self-adjoint  it  can  be  diagonalized  via  a 
unitary  transform.  If  we  use  the  Jacobi  method  [67]  to  obtain  the  eigenvectors  we  will 
gain  an  insight  on  the  effects  of  the  projection  on  these  eigenvectors.  We  can  use  the 
unitary  matrix  ,  v 

/cosi9,  sindi 


.   n  0   ,  (2.89) 

sm  Vi    —  cos  Vi 

to  diagonalize  (2.88),  and  obtain  its  eigenvectors  and  eigenvalues.  Once  i?,  is  determined 

the  eigenfunctions  are  obtained  from 

cos(i9,)  A,-  +  sin(i?,)  fii 


(2.90) 


sin(i?,)  A,  —  cos(i?,)  Hi 

and  i9,  is  determined  from 


^^2tf.- Wm+W£!M  (291) 


now 


{Xi\r\\i)^{fii\r\f^i)  =  Q  (2.92) 

and 

{^^\r\f^^)  =  {^^^\r\M)={l-Tf)^  (2.93) 
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then  matrix  (2.88)  is  of  the  form 


f    0  (i-r?)n 

v(i  -  T^)'    0  ; 


(2.94) 


which  has  eigenvalues  of 


±(i-r2)i 


(2.95) 


From  relations  (2.91)  and  (2.92)  we  see  that  tan2??j  =  oo,  or  t?i  =  for  all  pairs 
{i,  i').  This  leads  to  ... 


being  the  eigenfunctions  of  (2.88),  or,  the  same  eigenvalues  and  eigenfunctions  found 
previously  in  relations  (2.78-2.79)  and  (2.83-2.84). 

Harriman  solved  the  matrix  elements  of  (2.91)  for  the  projected  spin  density  operator^ 
and  obtained  i9  =  the  same  as  for  the  unprojected  density.  Thus  the  {<^f }  and  {<?i>f } 
are  unchanged  by  the  projection,  as  are  the  singly  occupied  {<f>'^}.  (The  interested  reader 
is  referred  to  References  62,  66  and  68  for  a  detailed  discussion  of  spin  projection  on 
a  UHF  wavef unction.) 

We  give  Harriman's  relations  for  the  projected  spin  density  matrix  eigenvalues  below. 
These  are  the  relations  used  in  this  work  to  explore  the  effects  of  spin  projection  on  the 

6.  Table  IV  of  Reference  66.  Harriman  solved  for  the  projected  eigenvalues,  and  eigenfunctions, 
of  p,  p+  and  p~  as  well  as  for  jo^. 


cosdi  sini?,- 
sin  t?i    —  cos 


(2.96) 
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calculated  aisn- 


(2.97) 


The  relations  (2.97)  are  for  the  special  case  of  l^l  =  5^  [66].  The  ^,  and  6'  are  paired 
eigenvalues  that  in  the  unprojected  state  sum  to  zero,  while  the  have  no  pairing 
properties  and  in  the  unprojected  state  are  equal  to  unity.  The  ujs{i)  term  is  defined  as 

g-l 

^s{i)  =  5](-l)^Cfc(5,m„f  -  1)  Afc(0  (2.98) 

A:=0 

where 

(V)  it 

Mi)  =  E  fl^^D  (2.99) 
and  the  Tj  are  those  of  equation  (2.36).  The  uj^  term  was  defined  in  relation  (2.66). 

It  was  not  our  intent  to  derive  the  relations  of  equation  (2.97),  to  do  so  would  fill 
half  a  Dissertation  with  unoriginal  work.  Harriman  and  others  [41,  51,  56,  57,  60-62, 
66,  68-73],  and  references  within,  have  worked  through  the  tedious  algebra  and  arrived 
at  (2.97).  We  have  implemented  (2.97)  for  the  expressed  purpose  of  investigating  the 
effect  of  spin  projection  on  the  a.^o  calculated  for  the  large  molecular  systems  we  are 
interested  in  describing. 
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The  AUHF  and  EUHF  Methods 

There  are  two  methods  related  to  PUHF  that  we  have  not  mentioned,  and,  except 

for  the  following  example,  will  not  pursue.  These  are  Annihilated  UHF  (AUHF)  [74, 

75]  and  Extended  UHF  (EUHF)  [70,  76].  In  Annihilated  UHF  one  or  more  of  the  higher 

multiplicity  components  are  removed  from  the  wavefunction  tfirough  the  application  of 

the  annihilator  Ag  [41].  As  in  PUHF  the  As  is/are  applied  after  the  SCF,  and  thus  the 

annihilated  wavefunction  is  no  longer  stationary  on  the  UHF  PE  surface.  In  practice  it  is 

usually  the  25"^  +  2  component  that  is  removed,  and  it  is  hoped  that  the  higher  multiplicity 

components  lie  too  high  in  energy  to  contribute  significantly  to  the  spin  contamination 

[43].  We  have  applied  the  Ag+i,  s  -  Sz,  to  the  INDO/l-UHF  wavefunction  describing 

the  Bacteriochlorophyll-a  (Bchla)  and  Bacteriochlorophyll-6  (Bchlfe)  cations  [77],  with 

the  following  results  for  (S^)  and  the  aN  of  the  four  pyrrole  nitrogens. 

Table  2-1:  The  affects  of  annihilation  on  {S")  and  calculated  a.so. 

Bchla+  Bchl6+ 


Experiment^ 

UHF 

AUHF 

1  Experiment 

UHF 

AUHF 

0.750 

3.970 

13.36 

0.750 

4.268 

14.79 

-0.78 

-2.24 

-4.68 

-0.82 

-1.80 

-8.75 

aN2 

-1.12 

-18.4 

-38.3 

-1.17 

-16.2 

-78.4 

aiVa 

-0.83 

-12.1 

-25.4 

-0.82 

-10.6 

-51.6 

aN, 

-1.03 

2.73 

5.52 

-1.06 

2.43 

11.5 

7.  From  Reference  77. 

8.  The  units  are  Gauss,  and  the  aN  are  for  ^*N. 
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From  Table  (2-1)  we  can  see  the  assumption  that  the  higher  multiplicities  contribute 
littie  to  tiie  total  wavefunction  is  wrong.  As  {S^)  increased  from  3.970  to  13.36  upon 
removal  of  the  |5|  =  |  component  of  tiie  INDO/l-UHF  wavefunction  for  Bchla+,  and 
from  4.268  to  14.79  for  Bchlfe"^.  Furtiier  we  see  tiiat  the  INDO/l-UHF  calculated  a.^o 
are  all  too  large,  and  more  than  double  upon  annihilation.  It  is  clear  that  the  components 
with  multiplicities  greater  than  {23^  +  2)  are  contributing  significantly  to  '^uhf  and  the 
computed  a,so  reflect  this.  (We  should  note  that  for  small  molecular  systems  annihilation 
does  improve  tiie  {S^)  as  well  as  tiie  calculated  energy  [78,  79].)  Then  from  Table  (2-1) 
AUHF  is  fatally  flawed,  at  least  when  viewed  from  the  perspective  of  wanting  a  method 
applicable  to  large  molecular  systems. 

The  Extended  UHF  method  mentioned  previously  is  basically  a  spin-projected  unre- 
stricted fvmction  where  the  orbitals  are  chosen  to  minimize  the  energy  after  projection. 
The  EUHF  wavefunction  is  stationary  with  respect  to  its  variational  parameters,  and  the 
spin  densities  calculated  with  it  tend  to  lie  between  those  determined  via  UHF  and  PUHF 
[70].  The  EUHF  methods  are  computationally  cumbersome,  and  as  stated  above  the  re- 
sultant spin  densities  are  no  better  the  those  of  PUHF,  at  least  in  the  sense  of  removing 
the  unwanted  spin  contamination.  For  these  reasons  the  EUHF  method  was  rejected  as  a 
possible  tool  for  exploring  spin  distributions  within  large  molecular  systems. 
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Calculated  Results  for  the  UHF  and  PUHF  Methods. 


In  Table  (D-8)  we  present  the  results  for  calculating  the  aH  and  aN  of  several 
more  test  structures  using  the  parameters  of  Appendix  D,  and  the  unpaired  spin  densities 
produced  via  the  methods  discussed  in  this  chapter.  In  Appendix  D  we  found  that  the 
INDO/1  Hamiltonian  yielded  consistently  worse  results,  compared  to  those  of  the  INDO/S 
Hamiltonian.^  For  this  reason  we  chose  to  limit  our  presentation  to  the  INDO/S  results. 
Along  with  the  additional  test  structures  we  also  present  results  for  Bchla  and  Bchl6 
cations  and  anions,  as  well  as  the  Bphea  anion.  These  are  the  types  of  compounds  we 
want  to  be  able  to  describe  properly.  The  structures  were  obtained  from  the  Brook  Haven 
National  Laboratories  crystallographic  data  bank,  through  a  collaboration  with  Dr.  Jack 
Fajer,  their  number  system  can  be  found  in  Figure  (15).  We  note  that  the  hydrogens 
of  the  three  crystal  structures  were  not  determined  experimentally  but  were  placed  in 
position  via  a  best  guess  algorithm.  Upon  inspection  of  the  Bchl  structures  we  found 
several  "strange"  methyl  groups  and  repositioned  the  hydrogens  associated  with  them  in  a 
manner  consistent  with  chemical  common  sense.  All  calculated  results  will  be  presented 
in  units  of  Gauss  (G). 


9.   See  Tables  (D-3)  through  (D-5). 
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Table  2-2:  Calculated  results  for  Bchla. 


Magnetic 

Calculated  a„o 

System 

Center^o 

Position^  ^ 

UHF 

PUHF 

Bchla+ 

2.11 

62 

N 

21 

-0.35 

-0.37 

-0.78 

N 

22 

-6.10 

-6.02 

-1.12 

IN 

ZD 

-2.88 

-2.86 

N 

24 

1.61 

1.58 

-1.03 

H 

5 

5.84 

6.82 

0.84 

H 

10 

7.30 

8.54 

0.46 

H 

132 

-3.23 

-2.07 

-0.58 

H 

20 

-1.55 

-1.74 

0.46 

Bchla- 

2.56 

55 

N 

21 

-2.05 

-2.03 

±0.42 

N 

22 

1.42 

1.43 

2.32 

N 

23 

0.38 

0.40 

±0.19 

N 

24 

1.76 

1.72 

2.09 

H 

5 

6.01 

7.16 

-3.43 

H 

10 

-7.19 

-8.42 

-2.46 

H 

132 

-2.98 

-1.64 

-2.22 

H 

20 

-2.12 

-2.47 

<±0.18 

10.  The  isotopes  are  "iV  and  ^H. 

11.  See  Figure  (15). 

12.  Reference  77. 


Table  2-3: 

Calculated  results  for  Bchl6  and  Bphea. 

Magnetic 

Calculated  Ui.^ 

System 

Center 

Position^  ^ 

UHF 

PUHF 

a 

Bchl6+ 

2.44  59 

N 

21 

-0.74 

-0.75 

-0.82 

N 

22 

-9.43 

-9.23 

-1.17 

IN 

Zj 

-2.24 

-2.24 

-0.82 

N 

24 

1.51 

1.49 

-  l.UO 

H 

5 

7.42 

8.28 

0.89 

H 

10 

8.64 

9.59 

0.43 

H 

132 

-2.25 

-1.76 

-0.43 

H 

20 

-1.03 

-1.10 

0.43 

Bphea- 

2.47  58 

N 

21 

-0.58 

-0.55 

-0.43 

N 

22 

7.15 

7.07 

2.56 

N 

23 

-0.78 

-0.75 

-0.21 

N 

24 

1.83 

1.71 

2.19 

CH3 

2' 

3.75 

2.71 

2.53 

H 

5 

-8.82 

-9.59 

-3.02 

H 

8 

-5.31 

-4.17 

-0.65 

H 

10 

-8.10 

-8.97 

-2.85* 

CH3 

12' 

4.67 

5.30 

2.96* 

H 

20 

5.80 

6.58 

-2.47 

N-H 

21' 

0.89 

0.84 

0.35 

N-H 

23' 

1.10 

1.07 

0.17 

13.  See  Figure  (15). 

14.  Reference  77. 
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The  results  presented  in  Tables  (2-2)  and  (2-3)  are  somewhat  disappointing.  Neither 
method  reproduces  the  experimental  observations  very  well.  All  four  cases  were  rather 
badly  spin  contaminated,  as  implied  by  the  value  of  (5^).  While  the  PUHF  method 
is  an  exact  eigenfunction  of  5^,  the  results  of  this  method  still  paralleled  those  of  the 
UHF  method.  The  only  apparent  affect  of  the  projection  was  to  diminish  the  \p'^\  by 
approximately  |,  as  is  expected  by  the  form  of  equations  (2.97).  Obviously,  after 
projection,  the  of  the  two  methods  stayed  proportional  to  one  another.  This  is 
surprising,  and  disappointing  as  we  had  hoped  that  the  projection  would  yield  values 
in  better  agreement  with  experiment,  having  removed  all  spin  contamination  from  a 
highly  contaminated  UHF  wavefunction.  Neither  the  signs  of  the  splittings  nor  the  their 
magnitudes  were  properly  described  by  either  method.  At  least  not  in  a  consistent  manner. 
We  should  also  note  that  the  *  values  are  considered  to  be  suspect  by  those  who  took 
the  measurement. 

It  is  obvious  from  these  results  that  neither  the  UHF  or  the  PUHF  methods  are  capable 
of  properly  describing  the  unpaired  spin  distribution  within  these  molecular  systems. 


CHAPTER  3 
THE  ROHF-CI  APPROACH 

Introduction 

In  this  chapter  we  shall  treat  the  open-shell  electrons  differently  at  the  SCF  level 
than  we  did  in  Chapter  2.  In  Chapter  2  we  used  Pople's  DODS  approach  and  obtained 
a  single  determinate  wave  function  which  was  an  eigenfunction  of  S'z,  as  it  must  be,  but 
not  an  eigenfunction  of  5^  In  this  Chapter  we  shall  use  an  SCF  procedure  that  produces 
a  wave  function  that  is  an  eigenfunction  of  both  these  operators.  Again  as  in  Chapter  2  it 
will  help  if  we  outiine  the  details  of  this  method,  the  Restricted  Open-Shell  Hartree-Fock 
(ROHF)  method,  as  it  has  been  developed  and  implemented  by  Zemer  and  Edwards  [80]. 

Restricted  Open-Shell  Hartree-Fock 

Restricted  Open-Shell  Hartree-Fock  was  first  proposed  by  Roothaan  in  1960  [81]. 
Since  that  time  numerous  methods  have  been  proposed  to  address  the  problem  of  open- 
shell  molecules  at  the  SCF  level  of  theory.  A  sampling  of  these  methods  can  be  found 
in  References  [81-92].  We  outline  the  Zemer-Edwards  method  below. 

The  total  wave  function  is  constructed  from  a  linear  combination  of  Slater  determi- 
nates. The  coefficients  of  the  Slater  determinates  are  chosen  to  represent  a  particular 
spin  state  [93],  to  guarantee  that  this  is  an  eigenfunction  of  5^.  The  coefficients  are 
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determined  by  spin  symmetry.  Fock-like  operators  are  constructed  for  each  shell, 
i.e.  closed-shell,  1-st  open-shell,  2-nd  open-shell  etc.  The  constraint 

(<^.|F''-F''|<^,)=0  (3.1) 

is  enforced  and  guarantees  the  matrix  of  Lagrange  multipliers  is  Hermitian,  where  Fi" 
is  the  Fock-like  operator  for  the  /^-th  shell  and  F"  for  the  i/-th  shell.  The  Hermiticity 
of  each  Fock-like  operator  guarantees  an  orthogonal  set  of  molecular  orbitals  for  each 
operator.  Though  the  various  sets  of  molecular  orbitals,  one  set  for  each  open-shell,  will 
not,  in  general,  be  orthogonal  to  one  another,  as  each  has  been  constructed  from  the 
diagonalization  of  its  own  operator  F**,  or 

i 

but  in  general 

(<^ri<^;)  =  (^ri<)7^o  (3.3) 

A  concatenated  set  of  orthonormal  molecular  orbitals  is  guaranteed  in  this  procedure  by 
projecting  out  the  eigenvectors  of  each  successive  Fock-like  operator  in  the  orthogonal 
complement  of  all  those  proceeding  it. 

The  spin  density,  equation  (1.13),  determined  from  an  ROHF  wave  function  is 
everywhere  positive,  a  consequence  of  the  construction  of  the  ROHF  wave  function. 
Although  negative  spin  densities  are  generally  small,  this  prediction  is  clearly  wrong. 
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Direct  experimental  observation  shows  the  sign  of  the  unpaired  spin  density  can  be 
negative  as  well  as  positive.  We  have  lost  some  of  the  flexibility  we  had  with  the  UHF 
wave  function,  a  price  we  paid  to  obtain  a  simultaneous  eigenfunction  of  the  Sz  and 
5^  operators  at  the  SCF  level.  A  single  spin-adapted  configuration  state  function  lacks 
the  flexibility  to  describe  differences  in  unpaired  spin  distributions.  We  would  have  to 
construct  a  total  wave  function  of  many  such  "spin-adapted  configurations"  to  have  the 
degree  of  flexibility  necessary  for  a  proper,  or  realistic  description.  The  simplest  way  to 
construct  such  a  total  wave  function  is  through  the  method  of  Configuration  Interaction 


The  Configuration  Interaction  method  is  a  straight  forward  application  of  the  Ritz  [94] 
method  of  linear  variations  as  applied  to  the  determination  of  electronic  wave  functions. 
Again  we  desire  to  solve  the  Schrodinger  equation,  equation  (1.14),  and  again  we  must 
approximate  its  solution.  In  the  CI  approach  the  approximate  wave  function  is 
expanded  as  a  linear  combination  of  "configuration  state  functions"  (CSF)  in  an  analogous 
manner  as  the  LC AO-MO  approach,  equation  (1.42),  for  the  Hartree-Fock  molecular 
orbitals. 


(CI). 


Configuration  Interaction 


N 


(3.4) 


s 
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Where  the  {$}  are  the  configuration  state  functions  and  the  expansion  coefficients,  Cs, 
are  the  parameters  which  are  varied  to  make  E[^aj,^]  stationary.  N  is  the  number  of 
expansion  functions.  We  examine, 

=  /  (*<.p.|*<.p.)  (3.5) 

in  which  H  in  equation  (3.5)  above  is  the  electronic  Hamiltonian,  equation  (1.20).  This 
leads  to  the  general  matrix  eigenvalue  equation 

HC  =  SCe  (3.6) 

where 

Hst  =  {^s\mi)  (3.7) 

and 

=  (3.8) 

The  linearly  independent  eigenvectors  Cp,  the  columns  of  C,  can  be  chosen  to  be 
orthonormal  with  respect  to  S, 

N 

=  S  ^tp-^^tc.,  =  (3.9) 

s,t 

leading  to  the  simple  eigenvalue  equation 

Hcp  =  c^E,  (3.10) 
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If  we  label  the  corresponding  eigenvalues  such  that 

Ei<E2<E3  <E^  '     ''  (3.11) 

each  Ej,  is  an  upper  bound  [95,  96]  to  the  corresponding  eigenvalue  of  Hei.  Furthermore 
as  additional  configuration  functions  are  added  to  the  expansion,  equation  (3.4),  each 
eigenvalue         of  the  {N  +  1)  term  expansion  satisfies 

<  ^(^^+1)  <  ^c^')  (3.12) 
p— 1  —   p      —  p 

and  as  {<t}s]  approaches  completeness  the  {E^}  approaches  the  exact  eigenvalues  of 
Hei  from  above  [95]. 

Correlation  Energy 

The  difference  between  the  "Hartree-Fock  limit"  energy,  which  is  the  limit  ap- 
proached by  restricted  self-consistent  field  calculations  as  the  basis  set  approaches  com- 
pleteness, and  the  exact  solution  of  the  electronic  Schrodinger  equation  has  been  called 
the  correlation  energy  [51].  It  reflects  the  fact  that  the  Hartree-Fock  Hamiltonian  contains 
the  average,  rather  than  the  instantaneous,  interelecd-on  potential,  and  thus  neglects  the 
correlation  between  the  motions  of  the  electrons. 

It  is  customary  today  to  break  down  the  correlation  energy  into  two  subcategories, 
dynamical  and  nondynamical.  Dynamical  correlation  is  the  concept  of  correlated  elec- 
tronic motion,  while  nondynamical  correlation  effects  are  due  to  near  degeneracies  and 
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rearrangement  of  electrons  within  partially  filled  shells.  The  CI  method  recovers  the 
correlation  energy,  both  dynamic  and  nondynamic,  at  least  in  the  limit  of  a  full  CI.  (We 
will  explain  the  various  CI  expansions  later.)  While  the  inclusion  of  correlation  into  an 
ab-initio  wave  function  is  necessary  for  the  accurate  description  of  observable  properties, 
it  is  a  less  than  good  idea  for  semiemperical  methods.  The  reason  this  is  so  is  straight 
forward;  semiemperical  methods  are  parameterized.  These  parameters  include  ionization 
potentials  and  electron  affinities  as  well  as  completely  empirical  parameters  such  as  P^,^, 
equation  (1.82),  which  are  varied  so  as  to  reproduce  some  benchmark.  If  the  parame- 
terization is  based  upon  experimental  values  then  correlation  is  included  to  an  unknown 
degree.  An  example  of  how  correlation  is  built  in  can  be  seen  in  the  Mataga-Nishimoto 
formulation  [97]  of  the  two-center  two-electron       integral,  equation  (1.80) 


where 

Jaa  =  IPa  +  EA,  (3.14) 

the  sum  of  the  ionization  potential  and  the  electron  affinity.  The  experimentally  measured 
IP  and  EA  are  already  "correlated";  any  and  all  such  effects  or  interactions  are  already 
embodied  in  these  measured  quantities.  To  now  introduce  correlation  into  the  INDO-SCF 
wave  function  is  to,  in  an  indeterminate  way,  "double-count"  the  correlation.  The  amount 
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of  correlation  that  is  built  in  depends  on  the  level  of  theory  that  is  parameterized  and 

what  parameters  are  chosen. 

Single  Excitation  Configuration  Interaction 

We  now  need  to  address  the  configuration  interaction  method  used  in  this  work. 
The  CI  basis,  {$},  is  constructed  from  a  set  of  spin  molecular  orbitals  obtained  from 
a  INDO-ROHF-SCF  calculation.  The  expansion  functions  $  are  linear  combinations  of 
Slater  determinates  whose  coefficients  are  chosen  such  that  the  CSF's  are  eigenfunctions 
of  botii  Sz  and  S\  where  (5^)  and  (S^)  of  {$}  are  tiie  same  as  tiiat  of  the  INDO- 
HF  reference  function.  Such  configuration  functions  are  said  to  be  "spin-adapted".  The 
method  used  to  construct  these  spin-adapted  functions  and  the  algebra  imposed  by  their 
form  in  the  evaluation  of  the  spin  density  will  be  addressed  in  detail  later.  We  now  need 
to  examine  the  expanded  wave  fimction  in  more  detail.  We  can  write  equation  (3.4)  in 
a  less  general  and  more  insightful  form  as 

occ  vir  occ  vir  occ  vir 

=  co$o  +  E  E     +  E  E      +  E  E  ^^k^^H  +  •  •  •  (^.is) 

i      1  i,j    a,b  i,j,k  a,b,c 

where  $o  in  the  first  term  of  equation  (3.15)  is  the  ROHF  reference.  The  set  {$^}  of 
the  second  term  is  comprised  of  all  spin-adapted  configurations  generated  by  removing 
one  electron  from  an  occupied  spin  molecular  orbital  ^,  and  placing  it  in  the  virtual,  or 
unoccupied,  spin  orbital  tpa-  These  are  known  collectively  as  the  single  excitations.  The 
is  generated  by  removing  two  electrons  from  the  occupied  molecular  orbitals,  or 
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occupied  space,  and  placing  them  in  the  virtual  space.  Such  spin-adapted  configuration 
functions  are  also  known  as  double  excitations  or  doubles.  The  same  pattern  holds  for 
the  triple  and  higher  excitations.  If  the  summation  of  the  various  excitation  levels  in 
equation  (3.15)  is  carried  to  its  limit,  then  the  calculation  so  performed  is  said  to  be  a  full 
CI  calculation.  The  number  of  configuration  functions  A'^  used  in  the  construction  of  a 
full  CI  wave  function,  for  a  given  number  of  electrons  n,  orbitals  M  and  total  electronic 
spin  angular  momentum  S  is  given  by  Weyl's  formula  [98] 

which  also  equals  the  dimension  of  the  full  CI  Hamiltonian  matrix.  Obviously  full  CI 
calculations  are,  from  a  computational  point  of  view,  intractable  for  all  but  the  smallest 
of  molecular  systems.  For  this  reason  the  CI  expansion,  equation  (3.15),  is  most  often 
truncated  to  include  at  most  only  the  first  two  or  three  terms.  Even  with  this  reduction 
in  the  size  of  the  problem  it  is  not  unusual  to  reduce  it  even  further  by  limiting  the 
size  of  the  "active  space".  That  is  to  limit  the  number  of  molecular  orbitals  from  which 
electrons  are  removed  from  or  placed  into,  usually  centered  around  the  highest  energy 
occupied  molecular  orbital  (HOMO)  and  the  lowest  energy  unoccupied  molecular  orbital 
(LUMO),  or  the  "HOMO-LUMO  gap".  For  the  work  in  tiiis  thesis  the  CI  expansion, 
equation  (3.15),  is  truncated  at  the  second  term  and  all  single  excitations,  are 
included  in  the  CI  wave  function,  and  it  is  referred  to  as  a  "singles  excitation  CI"  or 
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a  CIS  calculation. 

We  can  now  come  back  to  the  problem  of  the  "double-counting"  of  correlation.  If 
we  construct  a  simple  trial  wave  function  comprised  of  the  closed-sheU  SCF  reference, 
$0,  and  one  singly  excited  configuration  function,  or 

*c/s  =  (Co$o  +  Cf^n  (3-17) 

and  proceed  to  evaluate  {^cis\H^,\'^cis)  for  its  lowest  eigenvalue  and  associated  eigen- 
vector we  obtain  the  matrix  eigenvalue  problem 

It  is  obvious  that  any  mixing  of  the  diagonal  elements  is  through  the  off-diagonal  element 
($o|^e,|$n  =  ma)  +       [(^'H"^)  -  (^H^a)]  (3-19) 

r 

or  its  adjoint  The  m-th  element  of  the  Fock  operator,  again  in  the  basis  of  spin  molecular 
orbitals,  is  given  by 

mm  =  {i\h\a)  +  Yl^{ir\ar)  -  {ir\ra)]  (3.20) 

r 

Notice  that  the  right-hand-side  of  both  equation  (3.19)  and  equation  (3.20)  are  identical 
and  thus 

{^oA,m  =  {Hm)  (3.21) 
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Now  by  definition  solving  the  Hartree-Fock  eigenvalue  problem  requires  the  off-diagonal 
elements  of  the  Fock  matrix  to  vanish,  or  upon  solution 

■  '    •      ^  (^.im)  =  o  :;  ■  '  (3.22) 

and  then  so  must  equation  (3.19).  The  above  exercise  is  an  example  of  Bnllouin's 
Theorem  [99].  Which  is  stated  in  words  as;  singly  excited  Slater  determinates  will 
not  directly  interact  with  a  closed-shell  Hartree-Fock  reference  determinate  |$o)  through 
the  electronic  Hamiltonian,  or 

($o|^.,|$n=0  (3.23) 

as  well  as  its  adjoint.  Restricting  the  CI  expansion  to  the  {$"}  leaves  the  closed-shell 
Hartree-Fock  reference  uncorrelated,  and  so  the  problem  of  double-counting  correlation 
has  been  avoided.  This  approach  to  the  CI  method  has  proven  quite  successful  for  the 
calculation  of  the  low  energy  ultraviolet,  and  visible,  absorption  spectra  of  molecular 
systems  [100-103].  This  success  was  the  primary  reason  we  chose  this  technique  as  a 
starting  point  for  the  calculation  of  (jo^(r;^)). 

Unfortunately,  at  least  from  the  point  of  view  of  double-counting  correlation,  equation 
(3.23)  does  not  hold  for  an  open-shell  Hartree-Fock  reference.  Certain  single  excitations 
from  an  ROHF  reference  fimction  lead  to  spin-adapted  configuration  functions  which 
differ  from  the  reference  by  two  spin  molecular  orbitals  and  such  configuration  functions 
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interact  directly  with  |$o),  and  in  general 


($3|^.,|$^/)^0  (3.24) 

An  example  can  be  given  by  letting  the  wave  function  equation  (3.37)  interact  with  the 
wave  function  equation  (3.41)  through  the  electronic  Hamiltonian,  vis 

{T,\H^,\T,)  =  -^{im\ma)  (3.25) 
V6 

It  is  this  direct  mixing  of  the  reference  function  with  "certain"  singly  excited  functions, 
and  through  them  indirectly  with  the  other  single  excitations,  that  gives  the  CIS  approx- 
imation to  the  molecular  ground  state  the  flexibility  to  describe  differences  in  the  alpha 
and  beta  spin  distributions;  unlike  the  ROHF  reference  function. 

This  leads  to  the  question;  can  a  method  that  was  carefully  crafted  to  avoid  such 
interactions  as  equation  (3.25)  above  now  be  expected  to  well  describe  another  property, 
and  can  the  introduction  of  correlation  to  an  already  "somewhat"  correlated  reference 
function  be  done  in  such  a  way  as  to  obtain  consistent  results?  The  first  question  above 
is  answered  by  this  thesis  and  the  second  is  answered  in  Table  (3-1)  below. 
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Table  3-1:  How  the  correlation  of  the  INDO-ROHF  reference  affects  the  calculated  value  of  a.jo. 


% 

Correlation 

System 

aH  ExD^ 

A  flH 

CIS 

CISD^ 

^OH 

-26.4 

3.34 

8.31 

100.0 

-23.0 

2.49 

7.47 

98.4 

-20.2 

1.98 

4.23 

98.6 

2BH3- 

-15.1 

2.05 

5.50 

98.4 

-25.1 

3.53 

17.7 

100.0 

2NH3+ 

-25.9 

3.10 

9.99 

96.0 

2hcn- 

137.1 

4.14 

3.92 

93.5 

202H 

-9.8 

1.35 

1.38 

96.7 

In  Table  (3-1)  the  (aH  Exp)  column  lists  the  experimental  Hydrogen  aiso  values,  and  (A 
aH)  represents  the  difference  in  the  Hydrogen's  calculated  isotropic  hyperfine  splitting 
constants  when  calculated  from  a  CIS  and  a  CISD  wavefunction.  The  (%  Correlation) 
columns  lists  the  fraction  of  the  correlation  energy  recovered,  in  percent,  when  compared 
to  a  full  CI.  It  is  of  interest  to  note  the  average  (A  aH)  is  2.75  and  its  standard  deviation 
is  0.93;  recall  that  a,so  varies  almost  linearly  with  respect  to  p^,  equation  (1.10).  A 
comparison  of  columns  (aH  Exp)  and  (A  aH)  indicates  that  regardless  of  the  magnitude 
of  aiso  the  shift  is  a  rather  constant  3  Gauss.  Inspection  of  the  CIS  and  CISD  (% 
Correlation)  columns  shows  that  the  CIS  wave  function  has  recovered  a  small  fraction 


1.  All  aiso  values  are  reported  in  units  of  Gauss. 

2.  CISD  calculations  include  the  first  two  terms  of  equation  (3.15)  in  the  construction  of  the 
approximate  wave  function. 
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of  the  correlation  energy,  which  is  expected  due  to  the  generalized  Brillouin  theorem 
[104]  while  the  CISD  wave  function  has  recovered  almost  all  of  it.  Nevertheless,  the 
calculated  property  a,so  is  for  all  intents  and  purposes  the  same  for  either  wave  function. 
This  demonstrates  that  the  calculated  CIS  a.^o  values  increase  in  magnitude  by  a  small 
and  almost  constant  amount  when  the  wave  fimction  is  fully  or  almost  fully  correlated. 
This  small  "constant"  shift  indicates  that  a  CIS  wave  fimction  has  the  same  flexibility  to 
describe  the  spin  distribution  as  a  more  highly  correlated  wave  function  at  least  within  the 
INDO/S  model.  This  to  a  large  extent  relieves  our  concern  with  the  possible  problem  of 
correlating  an  already  correlated  reference,  and  the  possible  inconsistencies  that  may  have 
been  introduced  into  our  calculated  a.so-  We  chose  Hydrogen  for  the  above  endeavor 
for  the  following  two  reasons;  1)  there  is  no  core  orbital  spin  polarization  problem  when 
calculating  a,so  using  a  valence  orbital  only  theory  and,  2)  minimum  basis  set  theories 
have  problems  describing  the  "breathing"  of  the  atomic  orbitals  as  a  function  of  the  atomic 
charge.  By  taking  the  difference  in  the  computed  aiso  the  magnitude  of  the  contact  has 
to  a  very  large  extent  been  removed,  and  the  difference  represents  almost  purely  the 
difference  in  the  spin  density  as  computed  by  the  two  wave  functions.  This  leaves  (A 
aH)  a  good  gauge  of  the  effect  of  correlation  on  the  calculated  a,ao,  even  when  using  a 
minimum  basis  valence  orbital  only  theory. 
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Spin-Adapted  Configuration  Functions 

There  are  many  ways  to  construct  simultaneous  eigenfunctions^  of  Sz  and  5^  when 
starting  from  a  set  of  Hartree-Fock  spin  molecular  orbitals  {V*.  }-  Each  method  entails 
a  unique  procedure  for  the  evaluation  of  the  matrix  elements  of  the  CI  Hamiltonian,  a 
subject  we  shall  not  be  concerned  with,  as  well  as  for  any  integrals  necessary  to  evaluate 
the  expectation  values  of  observable  properties,  a  subject  that  shall  concern  us.  In  this 
work  two  Configuration  Interaction  procedures  were  used  in  the  evaluation  of  a,ao.  One 
is  quite  general  and  in  principle  can  treat  molecular  systems  of  any  multiplicity  and  of  any 
open-shell  structure;  it  is  also  quite  slow.  The  second  CIS  procedure  is  specific  to  a  single- 

— «  1 

reference,  singles-excitation,  calculation  on  molecular  systems  with  [Sz)  equal  to  |,  or 

doublet  molecular  systems;  it  is  also  quite  fast.  Since  at  the  time  of  their  implementation 

the  procedures  for  evaluating  (^(r^)}  for  either  CIS  method  were  unique,  we  shall 

present  them  in  some  detail  here. 

Rumer  Bonded  Functions 

A  generalized  extended  Rumer  diagram  methodology  [105-107]  has  been  presented 
by  Cooper  and  McWeeny  [108,  109].  Their  formulation  was  general  to  any  matrix 
element  between  bonded  functions  [106,  107]  and  included  formulae  for  spin-dependent 

operators.  We,  while  not  being  cognate  of  this  work,  also  worked  out  formulae  for  matrix 

3.  An  outstanding  reference  on  the  subject  of  "Spin  Eigenfimctions"  can  be  found  in  the  book 
by  Pauncz,  Reference  93. 
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elements  of  a  spin-dependent  operator,  built  upon  the  initial  work  of  Zemer  and  Manne 
[110].  Since  it  was  this  algebra  that  was  implemented  we  shall  develop  it  below  and  not 
that  of  Cooper  and  McWeeny,  through  as  one  would  expect  they  are  similar. 

The  valence-bond  type  spin  functions  grew  out  of  the  work  of  Heitler  and  Rumer 
[105-107]  and  of  Slater  [111],  and  are  commonly  called  Rumer  or  Slater  bonded 
functions.'*  The  configuration  functions  constructed  to  be  spin  eigenfunctions  via  the 
Rumer  method  are  linearly  independent  [93],  but  not  in  general  orthogonal.  The  Rumer 
functions  constructed  from  the  same  set  of  spatial  orbitals  will  not  be  orthogonal  among 
themselves  but  wiU  be  orthogonal  to  all  other  Rumer  functions,  assuming  an  orthonormal 
set  of  spatial  orbitals  form  the  initial  orbital  product.  The  spin  coupling  in  a  particular 
configuration  function  based  on  the  orbital  product 

M.---<l>n  (3.26) 
is  described  by  a  system  of  parentheses 

(^1  i<l>2  (^3  •  •  •  0„-.)^^„-i)^„)  (3.27) 

singlet  spin  coupling  pairs  of  orbitals  in  the  product.  For  a  given  number  of  electrons 
A'^  and  a  total  electronic  spin  angular  momentum  S,  there  will  be  y  -  5  pairs  of  spin 

4.  We  shall  use  the  notation  of  Boys,  Reeves  and  Sutcliffe  [112,  113]. 
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coupled  orbitals  for  the  case  of  5  =  nis.^  For  i  /  j 

{Mi)  =  M,{vMA^2)-^H^rWM  -  ^MccM]  (3.28) 

for  i  —  j 

(M,)  =  U.irMr.nc'MPM  (3.29) 

The  unpaired  orbitals  are  represented  by  "dangling"  left  parentheses,  again  for  i  /  ; 

{<j>icf>j){<t>>.  =  A<^.(rJ^,(r3)<^,(r3)-^[a(ti;,)/3(a;,)a(ii;3)  -  ^{u;,)a{u,)a{w,)]  (3.30) 

During  orbital  rearrangements  to  obtain  maximum  orbital  coincidence  between  two  spin- 
adapted  configuration  functions  each  orbital  keeps  its  associated  parentheses. 

Diagrammatically  the  Rumer  spin  functions  are  constructed  in  the  following  manner; 
the  spatial  orbital  labels  are  placed  at  the  vertices  of  a  polygon  with  one  extra  vertex 
added  as  the  "pole".  Starting  at  the  pole  and  proceeding  around  the  polygon  in  a  counter 
clockwise  manner  the  orbital  labels  will  be  in  ascending  order.  The  labels  are  then 
connected  by  arrows  into  singlet  coupled  pairs,  equation  (3.28),  with  the  arrow's  tail  at 
the  lower  index  and  its  head  at  the  higher  index.  As  many  such  couplings  as  necessary  are 
formed  to  obtain  the  desired  Sz  value.  All  uncoupled  orbitals  are  cormected  to  the  pole  by 


5.  For  a  given  N  electrons  and  S  electronic  spin  angular  momentum  there  will  be 
/iV'\     ?  N'  \ 

f{N',S)=  I,  j  -  spin  functions  [114]  for  a  given  orbital  product.  Where 

N'  equals  the  number  of  singly  occupied  orbitals  and  q'  equals  ^  -  S. 
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arrows.  Within  the  above  rules  all  such  possible  diagrams  are  constructed  and  any  with 
crossed  arrows  are  rejected.  Those  with  crossed  arrows  are  not  linearly  independent  [115] 
and  can  be  constructed  as  a  linear  combination  of  those  diagrams  whose  arrows  are  not 
crossed.  For  the  case  of  a  three  electron  three  orbital  S  =  ^  system  see  Figure  (2)  below. 
Equation  (3.31)  translates  the  diagrams  of  Figure  (2)  back  into  the  conventional  algebra 


1 


1 


1 


2  3 

Dl 


2  3 

D2 


Figure  2:  Rumer  bond  diagrams  for  the  case  of  three  electrons  in  three  orbitals  and  S  —  ^. 


Dl  =  {(f>i(l>i){<f>3  -  -^\(f>i(l)24>i\{oil3a  -  I3aa) 

D2  =  {M.){(f>,  =  -^\Ms<f>Ma/3a  -  ^aa)  (3.31) 


where  the  electron  labels  are  assumed  to  be  in  canonical  order  and  have  been  suppressed. 
By  inspection  it  is  obvious  the  Dl  +  Z)2  =  -D'i.  Armed  with  the  formalism  of  this 
section  we  are  ready  to  describe  the  evaluation  of  {p^{vn))  using  a  Rumer  CI  wave 
fimction. 
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Evaluation  of     for  a  Rumer  CI  Wavefunction 

The  evaluation  of  matrix  elements  betvi'een  Rumer  functions  for  the  spin  density 
operator,  p^,  is  achieved  as  follows.  Place  the  Rumer  functions  in  maximum  coincidence, 
then  connect  each  spin  molecular  orbital  in  the  first  function  to  its  counterpart  in  the 

second  with  a  vertical  bar.  vis 

{MM'l>.){<i>.  =  (1  2)  (5  (3  4) 

I   I     I    I   I  (3.32) 

{M2)i<i>M<l>.  =  {l  2)  (5  3)(4 
Then  horizontally  connecting  the  coupled  pairs  and  suppressing  the  indices  leaves  a 

"closed  loop"  and  a  "odd  chain",  vis 


Figure  3:  A  closed  loop  and  an  odd  chain. 

For  cases  where  S  is  greater  then  ^  "even  chains"  are  also  possible,  vis 


Figure  4:  Two  types  of  even  chains  which  arise  from  (12)  interacting  with  (1  (2  in  the  first  case 
and  (1  2)  (3  4)  interacting  with  (1  (2  3)  (4  in  the  second. 


and  in  general  so  are  multicycle  closed  loops,  vis 
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Figure  5:  A  multicycle  closed  loop,  arising  from  the  interaction  between  (1  2)  (3  4)  and 
(1  (2  3)  4). 

Once  the  matrix  element  between  a  pair  of  bonded  functions  has  been  reduced  to  the 
above  structures,  the  evaluation  of  (p^)  is  straight  forward.  The  operator  will  couple 
one  orbital  pair,  linked  by  a  vertical  bar,  at  a  time  and  {p^)  is  then  just  the  sum  of 
these  couplings. 

Now  from  the  properties  of  the  closed  loops  their  total  contribution  to  {p^)  sums  to 
zero.  This  is  a  direct  consequence  of 

{<l>,^,)  =  {<f>,ct>,)  (3.33) 

If  an  even  chain  is  formed  then  the  contribution  of  the  matrix  element  between  the  two 
Rumer  fimctions  is  zero.  There  will  always  be  spin  orthogonal  orbitals  coupled  in  such  a 
case,  and  no  one-electron  operator  leaving  the  primitive  spin  functions  unchanged,  such 
as  p^,  can  couple  even  chains.  This  leaves  only  the  odd  chains  contributing  to  {p^{rN)). 
Thus  only  the  density  contributions  of  the  orbitals  coupled  by  vertical  bars  in  the  odd 
chains  need  be  computed,  and  (^^(r;^))  as  evaluated  from  a  Rumer  CI  wave  function 
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is  written  as 


(3.34) 


I, J  Is 


or  in  the  atomic  orbital  basis 


{nr.))  =  E  "^^C-.r,,  Y.l'H  nAxMv  -  r,)|x.) 


(3.35) 


i,j  L  s 


where  the  summation  index  (5)  runs  over  all  the  vertically  linked  orbital  pairs,  9*  is  a 

link  parity  factor  and  alternates  as  H  h  •  •  •  +,  the  summation  indices  //  and  v  run 

over  all  atomic  orbitals,  the  coefficients  Cj  and  Cj  are  the  CI  expansion  coefficients,  see 
equation  (3.4),  for  the  /-th  and  J-th  Rumer  functions  respectively  and  the  indices  /  and 
J  run  over  all  such  bonded  functions.  The  atomic  orbital  density  P^^  is  determined  from 
the  5-th  vertically  linked  molecular  orbital  pair.  The  factor  Tjj 


accounts  for  the  product  of  the  normalization  factors  of  the  /-th  and  J-th  bonded 
functions  as  well  as  phase  changes  due  to  orbital  permutations  within  the  bonded  functions 
necessary  to  obtain  maximum  coincidence  between  them.  In  equation  (3.36)  (/)  is  the 
number  of  closed  loops,  (6)  is  the  number  of  bonded  pairs  (k)  is  the  number 

of  bonded  pairs  with  identical  orbitals  in  one  function  not  matched  by  a  bonded  pair 
of  identical  orbitals  in  the  other,  the  factors  (Ti  and  aj  are  the  number  of  permutations 
necessary  to  bring  the  unpaired  orbitals  back  to  their  original  order  in  bonded  functions 
/  and  J  respectively. 


T,,  =  {-2)'-\2)H-l) 


(3.36) 
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Fast  Doublet  CIS 


Whereas  the  Rumer  bonded  function  method  for  the  construction  of  spin-adapted 
configurations  was  very  general  with  respect  to  the  multiplicity,  number  of  open-shells 
and  degeneracies  within  the  open-shells,  the  following  method  for  their  construction  is 
very  specific,  and  entails  a  new  formalism  for  the  evaluation  of  (p^(r;v)).  For  lack  of 
a  better  name  we  shall  call  this  method  the  "Fast  Doublet  CIS"  (FDCIS)  method.  The 
FDCIS  algorithm  is  specific  to  the  following  case,  only  singles  excitations  off  a  single 
reference  determinate  of  spin  S  —  ^.  This  restriction  allows  for  concise  and  easily 
evaluated  matrix  elements  of  the  CI  Hamiltonian  [116],  again  this  topic  lies  outside  the 
scope  of  this  work.  We  shall  detail  below  the  formulae  for  the  evaluation  of  (p^(rjv)) 
from  a  FDCIS  wave  function. 

FDCIS  Spin-Adapted  Configurations 


Within  the  FDCIS  formalism  there  are  four  unique  molecular  orbital  configurations, 
see  Figure  (6),  these  lead  to  five  "types"  of  spin-adapted  configuration  functions.  We  shall 
label  these  generically  as  To  ■  ■  •  r4  and  shall  define  them  below.  In  an  attempt  to  achieve 
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in 


To  Ti  T2  T3&4 

Figure  6:  The  four  possible  orbital  configurations  for  the  FDCIS  formalism. 

compactness  of  form  the  following  convention  shall  be  followed;  a  Slater  determinant, 
l^i^i^^  ■  ■  ■  ^n\,  will  be  represented  as  \iij---n\.  The  orbital  indices  i,j,a  and  m  shall 
take  on  the  following  meaning,  i  and  j  will  represent  doubly  occupied  orbitals,  (a)  shall 
represent  a  virtual  orbital  and  m  shall  represent  the  open-shell  orbital  in  the  reference 
function.  Any  subscripts  on  the  orbital  indices  will  indicate  the  "type"  of  configuration 
function  from  which  it  came,  it  will  be  clear  below  what  this  means.  Now  let  us  write 
the  five  types  of  spin-adapted  configuration  functions 


To  =  '^o  =  \ll---m\ 


(3.37) 


T,  ='$7  =  |11  •• -mm- 


T,  =  '^l  =  \n-.-a\ 


T,  =       =  ■■■ia---m\  +  \ll---ai---m\) 


(3.38) 
(3.39) 
(3.40) 


83 

=       =  ■  ■  •  ia  •  •  •  m|  -  |11  •  •  •  ia  •  •  •  m|  +  |11  •  •  •     •  •  •  m|)  (3.41) 

where  the  superscripts  and  subscripts  on  $  indicates  the  single  excitation,  except  for  To 
which  is  the  ROHF-SCF  reference  configuration  function. 
Formulae  for  the  Matrix  Elements  of  Pf^cis 

Again  in  the  interest  of  compactness  we  shall  express  the  matrix  elements  of  the 
various  types  of  spin-adapted  configurations,  coupled  through  p^,  as  the  sums  of  the 
orbital  densities  they  contribute  to  p^^.   Let  us  start  with  the  diagonal  elements  and 


recognize  we  are  working  with  real  spatial  molecular  orbitals,  or  — 

{To\r\To)  =  \m){m\  (3.42) 

{TAr\Z)  =  \j){j\  (3.43) 

(r.l^^ir,)  =  \a){a\  (3.44) 

(r3|p^|r3)  =  \m){m\  (3.45) 

{T.\r\T.)  =  f         +  \a){a\  -  i|m)(m|)  (3.46) 

now  the  off-diagonal  elements  of  the  type  {Tj\p^\Tj)  where  Tj  ^  T'j 

{TAr\T[)  =  mj'\  (3.47) 

{T.\r\Tl)  =  \a){a'\  (3.48) 

(73  1^17;')=  0  (3.49) 
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{T.\r\T:)  =  i\iW\         ifa  =  a' 

(3.50) 

=  ||a)(a'|         if  i  -  i' 
now  the  off-diagonal  elements  of  the  type  {T^\p^\T„) 


Irri  1        1/77  \ 

=  biKm 

(3.51) 

{To\r%) 

(3.52) 

{To\r%) 

=  0 

(3.53) 

{To\r\T.) 

(3.54) 

=  0 

(3.55) 

{TAr\T.) 

if  ii  =  «3 

(3.56) 

{T,\r\T.) 

if  ii  =  «4 

(3.57) 

{nr\T.) 

=  -^V2){m\ 

(3.58) 

{TAr\T.) 

if  z'a       and  03  =  04 

(3.59) 

=  ^|a3)(a4|  if  ^3  =  ^4  and  03     04  (3.60) 

=  7^(|a)(a|  -  if     =     and  Ga  =  04 

This  leads  to 

{r{r.))  =      C*^^  E  Kt'^xMv  -  v,)\x.)  (3.61) 
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where  /  and  J  run  over  all  configuration  functions.  The  Ci  and  Cj  are  tiie  CI  expansion 
coefficients,  equation  (3.4),  tiie  indices  /x  and  v  run  over  all  atomic  basis  functions,  {x}. 
and  the  P'J  are  obtained  from  equations  (3.42)  through  (3.60). 

Calculated  Results  for  tiie  ROHF-CIS  Method 

As  in  Chapter  2  we  again  present  results  for  several  smaller  test  systems  similar 
to  those  used  in  determirung  the  optimum  valence  Slater  s-orbital  exponent,  in  Table 
(D-8).  As  was  the  case  in  Chapter  2,  the  INDO/1  Hamiltonian  yielded  consistentiy 
worse  results  than  did  the  INDO/S  Hamiltonian  when  compared  to  experiment.^  For  this 
reason  we  limit  our  presentation  of  results  to  those  of  the  INDO/S  Hamiltonian.  The 
Bacteriochlorophyll  structures  and  the  Bacteriopheophytin  structures  are  those  used  in 
Chapter  2,  and  again  all  calculated  results  are  presented  in  units  of  Gauss  (G). 


6.  See  Tables  (D-3)  through  (0-5). 
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Table  3-2:  Calculated  results  for  Bchla  and  Bphea. 


Magnetic 

Calculated 

System 

Center^ 

Position^ 

Bchla  + 

N 

21 

-0.09 

-0.78 

N 

22 

-0.47 

-1.12 

N 

23 

-0.52 

-0.83 

N 

24 

-0.08 

-1.03 

H 

5 

0.35 

0.84 

H 

10 

1.46 

0.46 

H 

132 

,  -1.36 

-0.58 

H 

20 

0.08 

0.46 

Bphea  - 

N 

21 

-0.62 

-0.43 

N 

22 

1.21 

2.56 

N 

23 

-0.09 

-0.21 

N 

24 

1.58 

2.19 

CH3 

2' 

2.47 

2.53 

H 

5 

-3.28 

-3.02 

H 

8 

-0.57 

-0.65 

H 

10 

-0.83 

-2.85* 

CH3 

12' 

0.51 

2.96* 

H 

20 

-3.67 

-2.47 

N-H 

21' 

0.15 

0.35 

N-H 

23' 

0.04 

0.17 

7.  The  isotopes  are       and  ^H. 

8.  See  Figure  (15) 

9.  Reference  77. 
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Table  3-3:  The  calculated  results  for  BcMft. 


Magnetic 

Calculated 

System 

Center 

Position 

a,so 

a  " 

Bchli  + 

N 

21 

-0.13 

-0.82 

N 

22 

-0.12 

-1.17 

N 

23 

-0.69 

-0.82 

N 

24 

-0.13 

-1.06 

H 

5 

0.24 

0.89 

H 

10 

0.64 

0.43 

H 

132 

-1.43 

-0.43 

H 

20 

0.23 

0.43 

Comparing  the  a„p  and  calculated  a„„  columns  of  Table  (3-2)  and  (3-3)  we  imme- 
diately see  that  the  ROHF-CIS  results  are  good.  The  computed  sign  of  the  hfc  are  all 
correct.  Unlike  the  UHF  and  PUHF  methods  the  computed  magnitude  of  the  hfc  are 
all  reasonable.  If  we  bear  in  mind  that  the  experimental  hfc  were  obtained  in  solution, 
and  the  structures  are  from  x-ray  crystallographic  methods,  then  these  results  do  indeed 
seem  reasonable.  We  should  also  note  that  the  *  values  are  considered  to  be  suspect  by 
those  who  took  the  measurement. 

While  it  required  obtaining  the  eigenvector  associated  with  the  lowest  eigenvalue  of 
a  matrix  of  order  26000,  these  results  seem  worth  the  effort. 


10.  See  Figure  (15) 

11.  Reference  77 


CHAPTER  4 
CONCLUSIONS 


Several  points  need  to  be  made  before  we  can  draw  our  final  conclusions.  Whether 
the  tools  we  use  to  investigate  nature  are  totally  empirical,  semiemperical  or  as  close  to 
first  principals  as  the  situation  will  allow,  there  is  only  one  real  measure  of  success;  and 
that  is,  will  our  model  describe  observed  reality?  This  leads  to  the  following  question, 
what  is  observed  reality?  Or  stated  in  a  less  philosophical  manner,  how  good  are  our 
measurements  of  nature?  The  qualification  of  this  question  will  shade  light  on  the 
"goodness"  of  our  model. 

Our  model  treats  a  single  molecular  system  as  the  sole  inhabitant  of  its  universe.  It 
interacts  with  nothing  but  itself  and  its  conformation  is  fixed.  Compare  this  to  the  world 
the  experimentalists'  molecular  systems  occupy;  most  experimental  observations  are  made 
in  the  solid  of  liquid  phases.  (All  molecular  systems  used  in  this  work  as  measures  to 
obtain  our  parameters,  or  as  yard  sticks  to  measure  the  parameters'  "goodness",  were 
themselves  measured  in  the  liquid  phase.)  In  the  few  cases  where  both  gas  phase  and 
condensed  phase  ESR  spectra  are  available  for  the  same  radical  species,  the  hfc  typically 
differ  by  a  few  percent,  but  can  differ  as  much  as  10+  %  for  the  extreme  cases. 

Computed  properties  such  as  spectroscopy,  dipole  moments,  transitions  moments, 
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etc.  as  well  as  ESR  hfcs  are  more  sensitive  to  structural  changes  in  the  molecular  system 
than  is  the  computed  energy  of  the  system.  The  hfc  are  particularly  sensitive  to  such 
changes  [117,  118].  Bartlett  et  al  [119]  in  a  work  to  be  published,  show  this  quite 
convincingly  for  triatomics,  at  several  levels  of  theory.  The  structures  we  used  in  the 
calibration  of  our  model  were  not  all  determined  experimentally,  many  were  obtained 
via  quantum  chemical  techniques.  The  conjugated  systems  were  optimized  via  an  INDO 
method,  where  practical  experience  has  shown  the  resultant  structure  to  be  as  good  as, 
or  better  than,  those  obtained  via  ab-initio  SCF  methods.  The  alkane  structures  were 
obtained  using  ab-initio  methods,  as  the  INDO  method  we  use  has  known  deficiencies 
when  applied  to  such  systems.  In  either  case  it  is  physical  intuition  and  low  level  quantum 
chemical  methods  by  which  the  structures  used  in  this  work  were  obtained.  Others  who 
work  with  smaller  systems  use  methods  which  are  capable  of  producing  more  precise 
structures,  usually  gas  phase  structures;  but  all  in  this  business  must  ask,  "is  the  structure 
we  are  using  the  same  as  Uiat  Uie  experimentalist  is  measuring"?  The  bacteriochlorophyll 
and  Bacteriopheophytin  structures  used  in  this  work  were  obtained  from  the  solid  state, 
and  the  hfc  ,  a„p,  are  from  solution  experiments.  The  "floppy"  substituents  of  the  Bchla, 
Bchl6  and  Bphea  are  prone  to  large  displacements  upon  going  from  the  liquid  to  the  solid 
state,  and  visa  versa.  As  noted  above  the  hfc  are  very  sensitive  to  changes  in  structure, 
so  we  chose  magnetic  nuclei  that  were  substituents  of  the  macrocycle  or  the  conjugated 
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rings  to  use  as  measures  in  our  calculations.  Here  we  assumed  that  these  rigid  regions 
would  be  the  least  perturbed  upon  phase  transitions. 

Molecular  vibrations  effect  the  measured  hfc  and  corrections  for  this  are  typically 
a  few  percent  of  the  measured  hfc  value  [118].  In  extreme  cases  where  large  atomic 
displacements  are  possible,  and  especially  when  such  displacements  move  the  measured 
magnetic  center  out  of  a  nodal  plane  into  direct  interaction  with  the  orbital  describing 
the  unpaired  electron,  such  as  the  case  of  the  methyl  radical,  the  vibrational  averaging 
correction  can  be  of  the  order  of  the  magnitude  of  the  observed  hfc  splittings.  The 
differing  slopes  of  Figure  (10)  may  be  partially  due  to  such  vibrations.  Vibrational 
modes  that  would  bring  the  hydrogen  lying  in  the  nodal  plane,  out  of  the  plane  and 
into  direct  interaction  with  the  tt  molecular  orbital  describing  the  unpaired  spin,  would 
increase  the  unpaired  spin  density  on  the  hydrogen.  Though  the  effect  would  tend  to 
reverse  the  sign  of  the  spin  density  at  the  hydrogen;  as  the  induced  spin  polarization 
at  the  hydrogen,  due  to  the  unpaired  electron  in  the  p^^  orbital  of  the  adjacent  carbon, 
would  have  a  negative  sign.  Direct  interaction  between  the  hydrogen  and  the  tt  electronic 
system  of  the  molecule  would  place  positive  spin  density  at  the  hydrogen. 

Multi-center  exchange,  which  is  missing  from  the  INDO  Hamiltonian,  would  of 
course,  also  effect  the  slopes  of  Figure  (10),  and  this  is  something  we  will  investigate 
in  the  future.  A  comparison  of  the  CIS  cr  and  tt  ^i,  of  Tables  (D-6)  and  (D-7)  indicates 
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that,  at  least,  for  the  systems  included  in  the  exponent  optimization  the  slopes  should 
be  very  close. 

Figures  (12-14)  are  plots  of  a,,„  vs  a,,p,  for  the  UHF,  PUHF  and  the  ROHF-CIS 
methods.  The  slopes  of  these  plots  are  determined  via  a  linear  least  squares  fit  and  should 
equal  1.0;  the  intercepts  should  equal  0.0.  The  UHF  and  ROHF-CIS  slopes  are  close 
to  1.0  (i.e.  1.01  and  1.04  respectively),  while  the  intercepts  are  respectively,  -0.99  and 
0.02.  While  the  slopes  for  both  methods  are  good,  the  intercept  of  the  UHF  method  in 
not  as  good  as  that  of  the  ROHF-CIS  method.  The  PUHF  slope  and  intercept  are  both 
bad,  compared  to  the  other  methods,  and  have  values  of  0.89  and  -1.39  respectively. 
Since  all  data  points  plotted  in  these  figures  are  for  systems  that  were  used  to  optimize 
^,  these  values  are  as  good  as  the  methods  are  capable  of.  From  above,  and  as  we  shall 
see  later  the  PUHF  method  is  the  least  successful  of  the  three,  followed  by  UHF,  and 
the  ROHF-CIS  method  will  prove  to  be  the  best. 

From  Figure  (10)  we  see  that  the  two-  and  three-center  terms  are  necessary.  If  we  do 
not  account  for  them  we  obtain  two  for  hydrogen.  One  for  the  a  systems  and  one  for 
the  TT  systems,  the  i„  1.48  while  1.20.  Accounting  for  the  two-  and  three-center 
terms  leads  to  =  1.248,  and  from  Figure  (14)  we  see  that  the  inclusion  of  the  two- 
and  three-center  terms  and  using  =  1.25  does  lead  to  good  results  for  hydrogen.  The 
aforementioned  exponent  being  that  of  the  ROHF-CIS  method,  since  both  the  UHF  and 
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PUHF  methods  proved  to  be  unsatisfactory,  see  Appendix  D. 

Other  researchers  have  shown  that  the  polarization  of  the  atomic  core  orbitals  can  be 
important,  as  well  as  electron  correlation  [118,  120].  Our  model  has  no  core  orbitals  and 
depends  on  the  core  effects  being  directly  proportional  to  the  valence  contribution  to  ai,„. 
This  has  been  shown  to  be  the  case  by  Chipman  [118]  for  the  second  period  elements 
when  the  structure  is  close  to  its  equilibrium  geometry. 

All  of  this  leads  us  to  ask  "just  how  close  to  the  experimental  values  can  we 
reasonably  expect  to  come"?  We  obviously  can  do  no  better  than  the  limits  placed 
on  a„p  by  interactions  with  the  environment.  Shifts  in  a„p  due  to  vibrations  are  not 
accounted  for  in  single  point  calculations.  Primarily  for  these  reasons  the  best  of  the 
high  level  quantum  chemical  methods  can  consistently  yield  results  of  no  better  than  ± 
20%  for  radical  molecular  systems.  It  would  then  seem  reasonable  to  assume  that  our 
model  should  on  average  target  ±  20%. 

The  results  of  Chapter  2  are  disappointing.  We  had  expected  the  removal  of  the  higher 
multiplicity  components  of  the  UHF  wavefunction  to  improve  the  calculated  hfc.  Instead, 
the  wavefunction  is  sufficiently  poor  even  after  spin  projection  that  no  improvement 
was  gained.  Neither  the  sign  nor  the  magnitude  of  the  calculated  hfc  splittings  can 
be  trusted  when  obtained  by  either  of  these  methods.  It  may  prove  interesting  to  take 
the  natural  orbitals  of  the  UHF  reduced  first-order  density  matrix,  which  are  unchanged 
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upon  projection,  and  use  them  and  their  projected  occupation  numbers  as  a  basis  for 
a  configuration  interaction  calculation,  or  a  perturbative  calculation,  to  obtain  the  spin 
density.  For  neither  wavefunction  is  useful  at  the  SCF  level  for  determining  a.^,,. 

The  results  of  the  ROHF-CIS  method  are  heartening.  The  signs  of  the  hfc  are 
correct  in  all  cases.  The  calculated  magnitude  of  the  splittings  are  reasonable.  We 
say  this  even  knowing  that  some  of  the  absolute  errors  are  large  (of  course  the  UHF 
or  PUHF  errors  were  far  larger).  The  Bacteriochlorophyll  structures  are  suspect,  and 
their  results  are  not  as  close  to  experimental  observations  as  one  would  have  liked.  The 
Bacteriopheophytin  structure  is  considered  as  good  as  one  can  obtain  for  such  molecular 
systems,  and  the  Bacteriopheophytin  calculated  results  are  noticeably  better  than  those 
of  Bacteriochlorophyll.  Several  points  need  to  be  made  or  reiterated  concerning  this 
point.  First,  all  hfc  are  quite  sensitive  to  structural  changes  in  the  molecular  system, 
and  there  is  little  chance  that  the  structures  we  use  in  these  calculations  were  identical 
to  those  measured  experimentally.  Secondly,  for  most  splittings,  and  for  all  negative 
splittings,  we  are  dealing  with  calculating  the  differences  in  small,  induced,  asymmetries 
in  spin  distributions.  When  calculating  small  hfc,  usually  for  centers  that  do  not  directly 
contact  the  unpaired  spin,  (e.g.  Those  magnetic  nuclei  that  lie  at  the  nodes  of  the  orbital 
describing  the  unpaired  spin.)  the  method  must  to  a  high  degree  be  capable  of  properly 
describing  both  the  alpha  and  beta  spin  distributions.  For  it  is  the  difference  in  these 
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two  distributions  that  determines  the  observed  hfc.  Seen  in  this  light,  this  method  is  at 
least  qualitatively  useful,  and  may  well  be  useful  as  an  aid  to  spectra  interpretation,  see 
the  results  for  Bphea-  in  Table  (3-2),  or  in  conformation  determination.  Finally  we  note 
that  measurements  on  most  large  systems  are  done  in  solution  or  in  condensed  phases. 
We  look  forward  to  modeling  this  effect  using  the  self-consistent  reaction  field  method 
[121],  which  experience  has  shown  to  have  major  effects  on  structure  [122,  123],  UV- 
Visible  spectroscopy  [124],  and  even  in  the  calculation  of  molecular  polarizabilities  and 
hyperpolarizabiUties  [125]. 

Finally,  it  may  prove  insightful  to  "classically"  give  these  molecules  a  temperature 
and  investigate  how  vibrations,  and  torsions  in  the  conjugated  system  affect  the  spin 
distribution.  Such  classical  mechanical  approaches  are  now  possible  to  apply  to  systems 
such  as  these. 


APPENDIX  A 
UNITARY  ORBITAL  ROTATIONS 

For  any  wavefunction  constructed  as  an  antisymmetrized  product  of  canonical 
Hartree-Fock  orbitals,  (in  general  any  orthogonal  functions)  the  expectation  values  de- 
termined from  such  a  wavefunction  are  invariant  to  unitary  rotations  of  the  canonical 
orbitals.  To  see  this  let  us  take  an  orthonormal  set  of  spatial  orbitals  {(f)}  and  with  them 
construct  a  square  matrix 

^M^)      •••  ^n(l) 

A=|     i       ••.       i     I  (A.1) 


of  order  n  x  n  such  that 


kM^)    •••  4>n{ri), 


^  =  -1=  det{K)  (A.2) 


Let  us  also  define  a  second  matrix  of  order  n  x  n  such  that 

=  (A.3) 
Such  a  matrix  is  said  to  be  unitary  and  when  used  in  transforming  the  {<j)}, 

ci>'i  =  J2<l>rUri  (A.4) 

r 

will  preserve  the  orthonormaUty  in  the  transformed  set.  Then  from  relation  (A.4)  it  is  clear 
that  the  matrix  which  corresponds  to  A  but  constructed  from  the  transformed  orbitals  is 

B  =  AU  (A.5) 
95 


96 

Now 

det{B)  =  det{AV)  =  det{A)det({J)  (A.6) 
which  is,  to  within  a  normalization  constant,  equal  to 

1^')  =  rfei(U)|*)  (A.7) 

Exploring  relation  (A.3)  leads  to 

det{U^U)  =  det{U^)det{U)  =  det{V)*  det(U)  =  \det({J)\^  =  det{l)  =  1  (A.8) 
or  in  general 

det{V)  =  e*^  (A.9) 

Then  ^'  differs  from  ^'  by  at  most  a  phase  factor,  and  if  U  is  a  real  matrix  this  phase 
factor  is  just  ±1.  It  is  a  postulate  of  Quantum  Mechanics  that  the  expectation  value  of 
any  observable  property  depends  on        thus  in  any  observable  way,  ^'  =  ^. 


APPENDIX  B 
THE  MIRROR  THEOREM 

Let  A  be  an  arbitrary  matrix  of  order  m  x  n,  and  let  B  be  an  arbitrary  matrix  of 
order  n  x  m  with  m  >  n.  The  products  AB  and  BA  are  quadratic  matrices  of  order 
m  X  m  and  n  x  n  respectively. 
Theorem: 

The  non-vanishing  eigenvalues  of  (AB)  and  (BA)  are  necessarily  the  same,  and 
the  classical  canonical  forms  associated  with  multiple  eigenvalues  are  identical 
[43,  53,  126]. 

The  mirror  theorem  as  stated  above  is  a  more  general  theorem  than  we  need.  We  choose 
to  demonstrate  the  theorem  for  the  special  case  of  product  matrices  that  are  self-adjoint, 
since  this  will  be  the  case  throughout  this  work.  A  more  rigorous  proof  of  the  general 
theorem  can  be  found  in  the  noted  references. 

If  S  is  Hermitian  (i.e.  =  S)  then  a  unitary  matrix  U  (i.e.  U+  =  U"^)  can  always 
be  found,  where 

U-^SU  =  A  (B.l) 

and  A  is  a  diagonal  matrix.  For  simplicity  we  assume  BA  to  be  Hermitian  and  have 
exactly  n  non-vanishing  eigenvalues.  Then  there  exists  a  matrix  V  of  order  n  x  n  that 
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brings  BA  to  diagonal  form 


V-^BAV  =  A  (B.2) 


Now  consider 


U  =  AVA-i  (B.3) 


where  we  choose  to  use  the  positive  square  root  of  the  A^.^  Then 

(AB)U  =  A(BAV)A-^  =  A(VA)  A-^  =  (AVA->)A 


(B.4) 


=  UA 

or 

U-^(AB)U  =  A  (B.5) 
Conversely  we  can  start  with  the  non- vanishing  eigenvalues  of  AB,  and  show  that,  if 

V-^ABV  =  A  (B.6) 

then  we  define 

U  =  BVA-^  (B.7) 

and 

(BA)U  =  B(ABV)A-^  =  B(VA)  A"^  =  (BVA-^)A  ^ 

(B.8) 

=  UA 


1.  This  forces  the  overlap  between  the  "paired"  orbitals,  Appendix  C,  to  be  positive. 


Equations  (B.2)  and  (B.5),  as  well  as  equations  (B.6)  and  (B.8),  are  eigenvalue  equations 
with  the  same  non-vanishing  eigenvalues,  A.  The  m-n  remaining  eigencolumns  of  AB 
could  be  repeat  eigencolumns  or  could  be  zeros.  Since  both  AB  and  BA  are  Hermitian 
and  both  have  a  rank  of  n,  each  then  have  n  orthogonal  eigencolumns,  and  the  m  -  n 
extra  columns  of  AB  are  then  zeros  and  so  have  zero  length. 


APPENDIX  C 
THE  PAIRING  THEOREM 

The  pairing  theorem  is  a  special  case  of  the  muTor  theorem,^  to  see  this  and  to 
demonstrate  the  theorem,  let  us  consider  two  sets  of  basis  fimctions, 

a^ia^aa,  ,a„,}  (C.l) 

b  =  {b,,b,  ,b„}  (C.2) 

of  order  m  and  n  respectively,  with  m  >  n.  Further  we  assume  the  sets  to  be  orthonormal. 
Then  one  has  the  overlap  integral  matrices 

S^,„  =  (a|b)        •  (C.3) 

St_  =  (b|a)  •  "  (C.4) 

Theorem: 

There  exists  two  unitary  transformations  U  and  V  having  the  property  that  if 

a'  =  aU 

(C.5) 

b'  =  bV 

then  the  overlap  matrices 

S'  =  (a'|b') 

(C.6) 

S't  =  (b'|a') 


1.  See  Appendix  B. 
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are  going  to  be  diagonal,  leading  to  a  natural  "pairing"  of  the  basis  functions. 
Consider  the  two  product  matrices  SS^  and  S^S,  which  are  Hermitian  and  of  order  mxm 
and  n  X  n  respectively.  (If  m  >  n  then  SS^  has  at  least  m  -  n  eigenvalues  equal  to 
zero.^)  Then  there  exists  a  unitary  matrix  V  such  that 

V-i(StS)V  =  u  (C.7) 

where  u  is  diagonal,  and  for  the  sake  of  simplicity,  assumed  to  be  of  order  n  x  n.  (i.e. 
No  vanishing  eigenvalues.)  Also  u„  >  0  since  S+S  is  positive  definite.  As  in  relation 
(B.3)  of  Appendix  B  we  define 

U  =  SVu~2  (C.8) 

where 

U^U  =  u-W^S+SVu-^  =  u"^  uu-i  =  (C.9) 

and 

U+SV  =  u- W+S^SV  =  u-^  u  =  (C.IO) 
If  for  convenience  we  define 

Tr=4  (C.ll) 

then 

  (a'rlb',)  =  TrSrs  (C.12) 

2.  See  Appendix  B. 
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as  well  as 

(bMa',)  =  TrSrs  (C.13) 

which  is  the  pairing  theorem. 

The  value  of  the  overlap  %  ranges  between  0  <  Tr  <  1.^  If  =  1  then  a'^  =  b'^ 
and  corresponds  to  a  doubly  occupied  orbital.  Those  orbital  pairs  having  Tr  <  1  are 
referred  to  as  the  "corresponding  orbitals"  of  Amos  and  Hall  or  the  "paired  orbitals"  of 
Lowdin.  Following  the  same  arguments  as  those  of  Appendix  B  leads  to  the  m  -  n 
extra  eigencolurans  of  SS^  associated  with  the  m  -  n  zero  eigenvalues,  all  necessarily 
being  of  zero  length. 


3.  As  long  as  the  positive  square  roots  of  the  Hk  are  chosen  in  relation  C.8. 


APPENDIX  D 
THE  EVALUATION  OF  \^pns{0)f 


In  Chapter  1  we  introduced  the  topic  we  wished  to  investigate  and  the  quantity 
we  intended  to  compute,  that  being  p^{r)  of  equation  (1.10).  In  Chapters  2  and  3  we 
investigated  methods  for  obtaining  a  "good"  P®;^  i.e.  good  in  the  sense  of  reproducing 
experimental  observations.  Here  we  shall  be  concerned  with  the  evaluation  of  the 
probability,  xjl(r)  x^lr).  of  equation  (1.89). 

Semiemperical  methods  in  general,  and  the  INDO/1  and  INDO/S  methods  we  inves- 
tigate in  particular,  never  explicidy  define  the  atomic  basis  {x}  [16];  rather  they  define 
Pock  and  other  energy  integrals.  This  leads  to  an  interesting  quandary;  upon  obtaining 
P'^,  or  P^,  how  does  one  go  about  obtaining  (O)  for  the  one-electron  operator  O, 


(where  O  operates  on  the  primed  basis  function,  and  1'  is  set  to  1  prior  to  integrating) 
when  one  does  not  have  a  basis  set?  Obviously  without  a  {x}  equation  (D.l)  is  difficult 
to  evaluate.  It  is  customary  in  semiemperical  theory  to  assume  a  basis  set  when  seeking 
to  evaluate  (D.l),  and  that  is  what  we  do  in  this  work.  We  select  as  our  basis,  a  set  of 

1.  See  equation  (2.63)  and  the  following  text. 


(D.l) 
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symmetrically  orthogonalized  [54]  Slater  type  atomic  orbitals  (STO)  [29], 

X„o  =  (Xl,X2.  Xm)  " 

(X..„|X,,„)  =  S  (D.2) 

or  LOwdin  type  orbitals  (LTO).  We  do  so  for  the  following  two  reasons,  1)  to  remain 
faithful  to  the  INDO  parameterization  [16]  and  2)  Slater  type  orbitals  behave  as  hydrogen- 
like atomic  orbitals  about  the  origin  [14],  i.e.  a  Slater  type  Is  orbital  is  identical  to  a 
hydrogen-like  Is  orbital.  All  hydrogen-like  s-type  orbitals  have  a  proper  cusp  at  the 
atomic  origin.  A  plot  of  the  nitrogen  Is  orbital  can  be  found  in  Figure  (7),  where  for 
Figures  (7),  (8)  and  (9),  r(x,  0, 0)  is  defined  in  the  following  way, 

r(x,y,z)  =  i/x2  +  y2  +  z2 

(D.3) 

r(x,0,0)  =  V^=  |x| 

The  unit  of  r  is  the  bohr  (oo)  which  is  equal  to  0.52918  x  10~^°  m,  the  Is  orbital  exponent 
used  in  Figure  (7)  is  taken  from  the  work  of  Clementi  [127],  and  has  the  value  of  6.6651. 

The  STO  radial  functions  are  nonorthogonal,  their  angular  parts  are  the  spherical 
harmonics  [14],  thus  for  orbitals  on  the  same  atomic  center  with  differing  angular 
momentum  quantum  numbers  /,  orthogonality  is  assured.  Though  those  with  the  same 
/  and  m  quantum  numbers  do  have  overlap.  The  normalized  set  of  Slater  orbitals  are 
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given  as 


x:*;,.(r)  = 


(2n)! 


2n+l 


r"-i  exp(-^r/ao)  Y,- 


(D.4) 


where  n,  /  and  m  have  their  usual  meaning,  (i.e.  The  principal  quantum  number,  n,  the 
orbital  angular  momentum  quantum  number,  /,  and  by  convention,  the  z  component  of 
the  orbital  angular  momentum,  m.)  The  normalized  spherical  harmonics  are  given  as 


Yr{0, 4>)  = 


"2/  +  1  {l-\m\)\ 


47r  (/+|m|)! 


P,'"*'  (cos^)  exp(— im^) 


(D.5) 


where 


(D.6) 


2'  /!  '  rfiy'+lmi 

and  It;  =  cos  0.  Only  5-type  atomic  orbitals,  I  =  m  —  0,  "contact"  the  nucleus  and  these 
have  no  angular  dependence.  From  equation  (D.4),  as  well  as  from  Figure  (7),  we  see 


that  the  l5  Slater  atomic  orbital  does  indeed  contact  the  nucleus,  with  a  value  of 

r  =  0  (D.7) 


where  Oo,  the  Bohr  radius,  is  defined  to  be  1  in  atomic  units. 


For  second  period  elements  the  value  of  xl',"  at  the  origin  is 


X',';{r)  =  N„  r  exp  (-^r)  =4>  0,      r  =  0 


(D.8) 


where 


2n+l  ■ 


iV„.= 


(20  

47r(2n)! 


(D.9) 
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this  is  true  for  all  Slater  3-type  orbitals  with  n  >  1,  see  equation  (D.4).  By  Schmidt 
orthogonalizing  the  25  to  the  Is  we  solve  two  problems,  1)  the  two  radial  functions  are 
now  orthogonal  and  2)  the  Slater  25  atomic  orbital  now  contacts  its  own  nucleus.  The 
orthogonal  2s'  is  given  as 

(1  -  ALJ^ 

where 

A,,..  =  (X:  J  X2,)  (D.ll) 

and  we  have  dropped  the  superscript  sto  from  the  {x}.  The  radial  functions  of  the 
hydrogen-like  atomic  orbitals  have  n  -  1  nodes.  It  is  obvious  from  (D.4)  that  the  STO 
radial  functions  have  no  nodes.  The  orthonormalization  has  introduced  the  single  node 
of  the  hydrogen-like  25  atomic  orbital  into  the  X2>'-.  Figures  (8)  and  (9)  are  plots  of 

2  2 

W2s-hyd{r)\  and  |X2.'('')|  ,  note  the  nodes  and  the  behavior  at  the  origin,  and  from 
Reference  14 

3 

V'2,-hvd(r)  =  r(  — )   (2-  — )exp(-Zr/2ao)  (D.12) 

where  Z  is  the  central  charge. 

In  Figure  (8)  Z  of  equation  (D.12)  has  a  value  of  5.0,  the  core  charge  of  a  Nitrogen 
atom,  assuming  perfect  shielding  by  the  two  I5  electrons.  The  Xz.'  valence  exponent, 
6„  is  that  of  nitrogen  and  has  the  value  of  1.95,  this  from  the  work  of  Zemer  et  al 
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[128,  129],  and  the  Is  exponent  is  the  same  as  that  used  in  Figure  (7),  i.e.  6.6651.  As 
expected  both  functions  are  very  similar  in  their  description  of  the  probability  density, 
especially  in  the  vicinity  of  the  origin.  This  is  the  region  we  need  to  well  describe  when 
calculating  the  Fermi  contact. 

The  above  diatribe  would  be  irrelevant  if  all  we  wanted  was  a  constant  to  multiply 
Tde/^egN^NP^,^,  by,  but  we  want  more  than  that.  We  chose  to  use  (D.7)  and  (D.IO) 
for  two  reasons,  1)  we  want  to  experiment  with  allowing  the  valence  s-orbital  contact  to 
vary  as  a  function  of  the  atomic  charge,  and  so  need  a  function  to  vary;  and  2)  we  want 
to  check  the  internal  consistency  of  our  model  by  using  the  observed  aiao,  the  {x}  and 
the  computed  Pf^^^  to  obtain  the  "observed"  valence  Slater  5-orbital  exponent  ^o6-  The 
closeness  of  (ob  to  ^p,  the  exponent  used  in  the  parameterization  of  the  INDO  model,  is 
a  measure  of  the  compatibility  of  the  INDO  generated  spin  density  and  the  "assumed" 
atomic  basis.  Tables  (D-1)  and  (D-2)  summarize  the  relationship  between  the  used 
in  the  parameterization  and  those  obtained  when  starting  from  the  experimentally 
observed  a.^o.  i.e.  The  ^o*  were  obtained  by  subtracting  the  computed  two-  and  three- 
center  contributions  to  a^„,,  from  the  observed  a.^o,  and  then  solving  equation  (1.10)  for 
^06-  The  inclusion  of  two-  and  three-center  terms  in  the  computed  a.so  has  the  effect  of 
slightly  improving  its  calculated  value  for  the  tt  radicals,  and  they  must  be  included  for 
the  a  radicals.  The  two-  and  three-center  contributions  are  large  for  a  radicals,  and  the 
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^06  of  the  cr  systems  is  far  too  large  if  only  the  one-center  terms  are  included.^  The  ^( 
of  both  a  and  tt  radicals  are  comparable  if  these  contributions  are  included.^ 


Table  D-1:  The     of  the  standard  Nitrogen  optimization  data  set,  where  ^,  =  1.95. 


System 

a^exp  (G) 

Charge 
(au) 

T  nwHin 

Charge 
(au) 

1,4,5,8-tetraazaanthracene  - 

±2.41 

1.935 

-0.434 

-0.348 

N,N-dihydropyrazine  + 

+  7.60 

1.972 

-0.222 

-0.020 

pyrazine  - 

±7.22 

1.985 

-0.532 

-0.467 

N,N-dihydrophenazine  + 

±  6.20 

1.919 

-0.254 

-0.071 

phenazine  - 

±5.14 

1.875 

-0.500 

-0.432 

N,N-dihydroquinoxaline  + 

±6.65 

1.921 

-0.236 

-0.048 

quinoxaline  - 

±5.64 

1.901 

-0.514 

-0.447 

The  Optimum  Slater  Valence  Exponent 


We  now  seek  the  optimum  valence  5-orbital  exponent  for  nitrogen  and  hydrogen.  We 
label  the  optimized  i  as  and  they  were  determined  in  the  following  manor.  A  standard 
set  of  molecular  systems  was  chosen  (All  experimental  information  can  be  found  in  the 
references  listed  in  Table  (D-10).)  and  their  structures  optimized  by  either  an  INDO/1- 
SCF  method  or  an  ab-initio  SCF  method  using  a  6-3 IG**  basis  set.'*  Then  the  aH  and 
aN  were  calculated  using  equations  (D.14)  and  (D.15)  respectively.  The      were  then 


2.  See  Table  (D-9). 

3.  See  the  CIS  values  of  Tables  (D-6). 

4.  See  Table  (D-10). 
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determined  such  that  the  error  between  the  calculated  ciiso,  Ocoio  ^nd  the  observed  ajso 
is  at  a  minimum.  The  function  minimized  was  the  root-mean-square  of  the  normalized 
errors,  i.e. 


/«)  = 


-y 


(D.13) 


N 

and  the  method  used  in  finding  the  minimum  was  that  of  Polak-Ribiere  [133,  134],  a 
gradient  based  method.  The  functions  used  for  aca/c(0  equation  (D.13)  are,  for  the 
valence  Is  and  the  2s  orbitals,  respectively, 

  "  V«..+6.)'-i92ff.a7 

where  C;^  is  a  constant  of  proportionality  specific  to  each  nucleus, 

:        C^  =  7jv(0.059621054)  (D.16) 

and  is  the  magnetogyric  ratio  [1]  of  nucleus  N.  The  numerical  constant  of  equation 
(D.16)  incorporates  the  constants  of  equation  (1.10),  as  well  as  conversion  factors  to  relate 
Jn,  expressed  in  Gaussian  units,  with  and  (  expressed  in  atomic  units,  to  yield  a  result 
in  Gauss.  The  nitrogen  of  (D.15)  is  again  that  of  Clementi  [127].  Equations  (D.14) 
and  (D.15)  were  obtained  by  the  evaluation  of  |V'i.(0)P  and  |V'2.'(0)|^  respectively  [135]. 
Tables  (D-3)  through  (D-5)  catalog  the  results  of  the  minimizations.  Before  we  review 
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the  results  we  need  to  comment  on  an  interesting  observation.  Figure  (10)  is  a  plot  of 
vs  Aiso,  where  A,so  is  defined  here  to  be  a^^p  minus  the  calculated  two-  and  three-center 
contributions.  It  is  obvious  from  Figure  (10)  that  the  cr  and  tt  hydrogens  are  treated 
somewhat  differently  by  the  INDO  Hamiltonian.  We  should  note  here  that  analogues 
plots  for  a  ab-initio  UHF  ST0-3G  and  3-2 IG**  basis,  yield  similar  results.  Part  of  the 
difference  is  due  to  the  two-  and  three-center  contributions.  The  curve  labeled  3CNT  is 
a  plot  of  a^^p  vs  /9f,  where  the  two-  and  three-center  terms  are  accounted  for;  the  curve 
labeled  a  does  not  have  these  terms  taken  into  account,  and  the  curve  labeled  tt  is  a  plot 
of  the  TT  radical  a^^^  vs  Aiso-  From  Table  (D-6)  we  see  that  the  CIS  hydrogen  a  and  tt 
are  quite  similar,  and  while  their  plotted  curves  still  differ,  the  (op  for  the  a  compounds 
has  changed  from  1.454  to  1.250  upon  accounting  for  the  two-  and  three-center  terms. 
Sigma  systems  interact  directly,  while  for  tt  systems  exchange  is  important.  Multicenter 
exchange  is  missing  in  the  INDO  Hamiltonian,  and  this  can  be  expected  to  have  an  effect 
on  the  c7  and  TT  slopes.  Figure  (14)  clearly  shows  that  once  die  two-  and  three-center 
terms  are  accounted  for,  and  is  used,  the  apparent  difference  between  the  cr  and  tt 
systems  is  no  longer  there.  We  will  address  this  further  in  the  conclusions.  The  results  of 
Tables  (D-6)  and  (D-7)  are  obvious.  For  hydrogen  and  nitrogen  the  INDO/S  Hamiltonian 
is  superior  to  the  INDO/1  Hamiltonian.  Within  the  INDO/S  model,  ROHF-CIS  gives  the 
best  results.  For  these  reasons  we  will  present  the  INDO/S  results  in  Chapters  (2)  and 
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(3),  and  refrain  from  presenting  the  INDO/1  results. 

It  is  of  interest  to  note  that  the  hydrogen  of  Table  (D-6)  approach  the  same  value 
for  (7  and  tt  radicals  for  the  CIS  and  PUHF  wavefunctions,  but  do  not  for  the  UHF 
wavefunction.  It  would  seem  the  UHF  wavefunction  proportionately  over  estimates  the 
spin  density  on  the  hydrogens  lying  in  the  nodal  plane  of  tt  radicals,  while  PUHF  and 
CIS  seem  to  get  it  proportionately  right. 
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Table  D-3:  Minimization  results  for  tt -hydrogen. 


Optimized  Valence  5-orbital  Exponent 


UHF 

PUHF 

CIS 

INDO/S 

(Error  %) 

40.2 

40.2 

6.14 

Max  Err  % 

-98.7 

-98.0 

-11.2 

INDO/1 

(Error  %) 

42.8 

42.8 

18.0 

Max  Err  % 

-96.5 

-96.4 

-32.5 

Table  D-4:  Minimization  results  for  cr-hydrogen. 

Optimized  Valence  5-orbital  Exponent 

UHF 

PUHF 

CIS 

INDO/S 

(Error  %) 

6.89 

7.55 

6.49 

Max  Err  % 

-13.0 

-13.1 

-12.4 

INDO/1 

(Error  %) 

26.7 

14.2 

18.0 

Max  Err  % 

-85.8 

-33.0 

-32.5 

Table  D-5:  Minimization  results  for  tt 

-nitrogen. 

Optimized  Valence  5-orbital  Exponent 

UHF 

PUHF 

CIS 

INDO/S 

(Error  %) 

8.92 

8.91 

4.07 

Max  Err  % 

15.3 

15.1 

-9.74 

INDO/1 

(Error  %) 

17.9 

8.91 

6.86 

Max  Err  % 

35.9 

15.1 

12.8 
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Table  D-6:  The  optimized  hydrogen  exponents. 


Optimized  Valence  5-orbital  Exponent 
UHF  PUHF 

CIS 

a  one-center 

top 

1.406  1.454 

1.424 

a  one-,  two-  and  three-center 

1.252  1.295 

1.274 

a  and  tt  one-,  two-,  and  three-center 

0.961  1.389 

1.248 

Table  D-7:  The  optimized  nitrogen  exponents. 

Optimized  Valence  5-orbital  Exponent 
UHF  PUHF 

CIS 

1.792  2.503 

1.956 

The  Charge  Dependence  of  i 


We  had  initially  hoped  to  find  a  relationship  between  the  computed  charge  of  the 
atomic  center,  either  the  Mulliken  charge  [130,  131]  or  the  Lowdin  charge  [132],  and  the 
"necessary"  |V'(0)p;  i.e.  the  necessary  to  obtain  a„p.  Such  a  relationship  would  have 
allowed  the  valence  5-orbital  to  expand  or  contract,  "breath",  as  a  function  of  the  charge 
associated  with  each  magnetic  center.  Inspection  of  Tables  (D-1)  and  (D-2)  clearly  shows 
that  there  is  no  such  relationship,  and  we  abandon  this  line  of  investigation  here. 
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Table  D-8:  Computed  UHF,  PUHF  and  ROHF-CIS  results  for  the  optimum  Slater  valence 
exponent 


Magnetic 

Molecular  System 

Center' 

T  TT  TT" 

UHF 

TtT  TT  TT" 

PUHF 

CIS 

phenazine  - 

0.847 

H 

-3.66 

-1.70 

-1.15 

±  1.93 

-0.45 

-0.21 

-0.70 

±  1.61 

N 

6.27 

6.11 

5.10 

±  5.14 

N,N-dihydrophenazine  + 

0.783 

H 

-5.38 

-5.44 

-7.38 

±  6.79 

-1.12 

-1.19 

-1.28 

+  1.77 

0.44 

0.30 

-0.23 

±  0.59 

N 

6.36 

6.32 

6.10 

+  6.20 

quinoxaline  - 

0.822 

H 

-1.85 

-2.03 

-3.24 

+  3.32 

-1.39 

-1.46 

-1.75 

+  2.32 

-0.16 

-0.15 

-0.74 

±  1.00 

N 

6.26 

6.08 

5.08 

±  5.64 

NJ>J-dihyclroqumoxalme  + 

0.767 

H 

-5.90 

-5.94 

-8.14 

+  7.17 

-2.56 

-2.61 

-4.67 

±  3.99 

-0.88 

-0.92 

-1.01 

±  1.38 

0.27 

0.34 

-0.26 

±0.78 

N 

6.41 

5.49 

6.52 

±6.65 

oxacyclobutane  + 

0.751 

H 

66.5 

75.5 

63.7 

63.0 

12.2 

14.2 

11.6 

11.0 

tetrahydrofuran  + 

0.751 

H 

86.0 

98.3 

82.0 

88.5 

44.2 

50.7 

42.1 

39.8 

7.  The  isotopes  are  "iV  and  ^H. 
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Table  D-9:  The  effects  of  including  the  two-  and  three-center  contributions  to  a 


Number  of  Included  Terms 

Magnetic 

1  Term 

3  Terms 

System 

Center 

flcaic  (G) 

acalc  (G) 

Oexp  (G) 

NO2 

N 

56.2 

55.1 

52.5 

PH2 

H 

-10.8 

-13.5 

-17.9 

CeHg  - 

H 

-4.95 

-4.89 

-3.75 

NH3  + 

N 

20.2 

21.1 

19.5 

H 

-21.9 

-25.5 

-25.9 

C4H8O 

H 

55.7 

82.0 

88.5 

H 

28.6 

42.1 

39.8 

CsHeO 

H 

44.7 

63.2 

63.0 

H 

8.37 

11.61 

11.0 
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Table  D-10:  The  Optimization  Data  Set. 


Chemical  Name 

Formula 

Charge 

Reference 

Optimized 

oxacyclobutane 

CsHeO 

+1 

136 

6-31G** 

tetrahydrofuran 

C4H8O 

+1 

137 

6-31G** 

anthracene 

C14H10 

-1 

138 

INDO/1 

naphthalene 

CioHg 

-1 

139 

INDO/1 

anthracene 

C14H10 

+1 

138 

INDO/1 

phenazine 

C12H8N2 

-1 

140 

INDO/1 

quinoxaline 

CgHeNa 

-1 

141 

INDO/1 

N,N-dihydropyrazine 

C4H6N2 

+1 

140 

INDO/1 

1 ,4,5,8-tetraazaanthracene 

C10H6N4 

-1 

141 

INDO/1 

pyrazine 

C4H4N2 

-1 

141 

INDO/1 

nitrogen  dioxide 

NO2 

0 

142 

6-31G** 

nitric  oxide 

NO 

0 

142 

6-31G** 

benzene 

CeHg 

-1 

142 

INDO/1 

ammonia 

NH3 

+1 

142 

6-3 IG** 

N,N-dihydrophenazine 

C12H10N2 

+1 

143 

INDO/1 

N,N-dihydroquinoxaline 

C8H8N2 

+1 

143 

INDO/1 

The  Results  of  Optimizing  the  Valence  Exponent 


In  Table  (D-8)  we  summarize  the  results  of  our  optimization  of  the  valence  5-orbital 
exponents.  It  is  clear  that  the  smaller  the  hfc,  the  greater  the  absolute  error.  Overall 
though,  the  results  are  quite  acceptable.  Of  course  the  real  test  will  come  when  we  apply 
the  results  of  this  appendix  to  the  Bacteriochlorophyll  and  Bacteriopheophytin  radical 
ions.  The  numbering  for  these  structures  can  be  found  in  Figure  (15). 

From  Table  (D-8)  we  see  that  while  the  magnitude  of  the  splittings  are  "reasonable" 
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for  all  these  methods,  the  signs  are  suspect  All  7r-radical  hydrogens  are  expected  to 
have  negative  splittings,  and  as  we  can  see  only  the  ROHF-CIS  method  predicts  this 
for  all  such  hydrogens.  As  we  shall  see  later  the  ROHF-CIS  method  always  correctly 
predicts  the  sign  right 


Figure  8:  The  Hydrogen-like  2s  probability  density. 
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Figure  9:  The  Slater  type  2s'  probability  density. 
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Figure  10:  A  plot  of  the  spectroscopic  QS  P'     vs  the  experimental  A.v 
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Rgurc  11:  The  x  region  of  Figure  (10)  expanded. 
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Figure  12:  A  plot  of  UHF  a„„  vs  a„,.  using  ^..=  0.961. 
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Figure  13:  A  plot  of  PUHF  a„,<,  vs  a„„  using     =  1.389. 
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Figure  14:  A  plot  of  a„,,  vs  a„,,  using  ^,,=1.248. 
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